ApJ, 750, 43; submitted: May 25th, 2011; accepted: February 25, 2012 
Preprint typeset using LAT^X style cmulateapj v. 5/2/11 



(N 

o 

(N 
u 

< 



Of 

6 

o3 



NEPTUNE'S WILD DAYS: CONSTRAINTS FROM THE ECCENTRICITY DISTRIBUTION OF THE 

CLASSICAL KUIPER BELT 

Rebekah I. Dawson 1 & Ruth Murray-Clay 

Harvard-Smithsonian Center for Astrophysics 
60 Garden St, MS-51, Cambridge, MA 02138 
ApJ, 750, 43; submitted: May 25th, 2011; accepted: February 25, 2012 

ABSTRACT 

Neptune's dynamical history shaped the current orbits of Kuiper Belt objects (KBOs), leaving clues 
to the planet's orbital evolution. In the "classical" region, a population of dynamically "hot" high- 
inclination KBOs overlies a flat "cold" population with distinct physical properties. Simulations of 
qualitatively different histories for Neptune including smooth migration on a circular orbit or scattering 
by other planets to a high eccentricity have not simultaneously produced both populations. We explore 
a general Kuiper Belt assembly model that forms hot classical KBOs interior to Neptune and delivers 
them to the classical region, where the cold population forms in situ. First, we present evidence that 
the cold population is confined to eccentricities well below the limit dictated by long-term survival. 
Therefore Neptune must deliver hot KBOs into the long-term survival region without excessively 
exciting the eccentricities of the cold population. Imposing this constraint, we explore the parameter 
space of Neptune's eccentricity and eccentricity damping, migration, and apsidal precession. We rule 
out much of parameter space, except where Neptune is scattered to a moderately eccentric orbit 
(e>0.15) and subsequently migrates a distance Aon = 1 — 6 AU. Neptune's moderate eccentricity 
must either damp quickly or be accompanied by fast apsidal precession. We find that Neptune's 
high eccentricity alone does not generate a chaotic sea in the classical region. Chaos can result 
from Neptune's interactions with Uranus, exciting the cold KBOs and placing additional constraints. 
Finally, we discuss how to interpret our constraints in the context of the full, complex dynamical 
history of the solar system. 

Subject headings: Kuiper Belt, planets and satellites: Neptune, solar system: general 
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1. INTRODUCTION 

Neptune, with its nearly circular and equatorial or- 
bit, may seem straight-laced compared to the oblique, 
hot, eccentric, and resonant planets in the extra-solar 
menagerie. B ut the highly inclined and eccentric or- 
bits of Pluto (|Malhotral 119931 [1995) and subsequently 
discovered Kuiper Belt objects (KBOs) imply that Nep- 
tune may have experienced its own "wild days" in the 
early solar system. During these wild days, Neptune 
sculpted the KBOs into four main dynamical classes: 
objects in mean motion orbital resonance with Nep- 
tune (the "resonant" population), objects that are cur- 
rently scattering off Neptune (the "scattering" popula- 
tion), and two populations of "classical" objects that 
are currently decoupled from Neptune. One population 
of classical objects is dynamically "cold," on nearly cir- 
cular orbits at low inclinations, and the other classical 
population is dynamically "hot" with a range of eccen- 
tricities and inclinations. The cold classicals have dis- 
tinct physi cal properties from the hot classicals, includ- 
ing colors dTegler fe Roman ishin 20 00t iThommes et al.l 
2001 iPeixmho et al.ll200 8l). sizes (ILevison fe Sternl l2001l: 
Fraser et al.ll2010D. albedos (iBrucker et alj|2009|). and bi- 
nary fraction ([Stephens fc Nol]||2006HNoll et al.ll2008| ). A 
major problem in understanding the formation of the so- 
lar system is that, as we will review below, no model 
of Neptune's dynamical history adequately produces the 
superposition of hot and cold classicals or accounts for 
the difference in their physical properties. 

1 rdawson@cfa.harvard.edu 



Two types of dynamical sculpting models have been 
developed to explain, in particular, the pop ulation of 
resonant KBOs. Extensive migratio n models ([Malhotral 
119931 119951 : lHahn fc Marhotrall2005f ) propose that Nep- 
tune migrated outward by 7-10 AU on a nearly circular 
orbit from its location of formation, capturing objects 
into resonance as its resonance locations slowly swept 
through the Kuiper Belt. This type of model generates 
the resonant and scattering objects and the cold popu- 
lation (unexcited objects that, in this model, formed in 
situ) but not the hot population. It also does not match 
the observed inclination distri bution within the r eso- 
nances. Chaotic capture mod els (ILevison etaLf2008), in- 
spired by the Nice model (sec Morbi delli et al.ll2008l and 
references therein), propose that Neptune was scattered 
onto a highly eccentr ic orbit by other plan ets during a 
period of instability ()Thommes et al.l 119991 ) . Neptune's 
high eccentricity created a chaotic zone in what is now 
the classical region, and some objects were caught in res- 
onances when Neptune's eccentricity damped. Chaotic 
capture models produce a hot population: objects that 
formed in the inner disk, subsequently were scattered 
by Neptune into the classical region, and then decou- 
pled when Neptune's eccentricity damped. These mod- 
els also produce a resonant population and scattering 
population. Although some of the objects delivered into 
the classical region end up on low-eccentricity orbits, we 
point out that a cold population confined to low eccen- 
tricities is not produced. Other variations of the Nice 
model fe.g. lMorbidelli et al.ll2008[ ) include an in situ pop- 
ulation of cold objects, but, over the course of Neptune's 
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evolution, these objects become excited to higher eccen- 
tricities. K. Batygin (2010, private communication^) has 
suggested that fast apsidal precession of Neptune's orbit 
could prevent Neptune from disrupting the cold classi- 
cal during its proposed high-eccentricity period, but it 
remains to be explored under what circumstances this 
mechanism would work and how it would affect the hot 
classicals. 

With neither the extensive migration models nor 
chaotic capture models producing both the hot and cold 
classicals, the qualitative picture of what happened in the 
early solar system, including the roles of planet-planet 
scattering and planetary migration, remains up for de- 
bate. It remains a question whether Neptune migrated 
outward by many AU on a nearly circular orbit, was 
launched onto an eccentric orbit near its current loca- 
tion, or none of the above. Pinning down Neptune's 
dynamical history, which should be possible given the 
constraints from over 500 KBOs with well characterized 
orbits, would reveal much about the history of our own 
solar system and about the processes of scattering and 
migration that shape the architecture of many planetary 
systems. 

Previous models attempted to produce all four dynam- 
ical classes of KBOs with 0.5-4 Gyr simulations that in- 
cluded all four giant planets and thousands of massles s 
KBOs (e.g. lHahn fc Malhotral[2005l: ILevison et al.H2008D . 
It has not been computationally feasible to fully explore 
parameter space with such extensive simulations. Thus it 
is unclear whether the dynamical history described by a 
particular model (1) has trouble producing both the cold 
and hot populations because there is a qualitative prob- 
lem with the scenario or, alternatively, because the pa- 
rameters need to be slightly adjusted; and (2) is unique, 
or whether another, qualitatively different dynamical his- 
tory would match the observations just as well. 

Inspired by previous models, we explore a general- 
ization in which Neptune undergoes all potential com- 
binations of high eccentricity, migration, and/or apsi- 
dal precession: "Neptune's wild days." In this gener- 
alization, the cold objects form in situ where we ob- 
serve them today and the hot objects are delivered from 
the inner disk and superimposed on the cold objects. 
"Two-origin" models superimposing a hot classical pop- 
ulation from the inner disk on top of a cold popula- 
tion f irmed in situ (e.g. | Levison fc Sternl l200il IGomesI 
120031 : iMorbidelli et al.ll2008f ) have the advantage of ex- 
plaining the different physical properties of the hot and 
cold classicals that were discussed above. The differ- 
ent colors, sizes, and albedos of the two populations are 
accounted for by their formation in different regions of 
the solar system's proto-planteary disk under different 
conditions. For instance, chemical differences may re- 
sult in different colo rs for objects form ed in the inner 
versus the outer disk (jBrown et al.ll20lil ). The cold clas- 

2 After the submission of this manuscript, Batygin ct al. (2011) 
presented a model in which Neptune underwent a period of high 
eccentricity and, due to its fast apsidal precession, could avoid dis- 
rupting the cold classicals. Because this paper appeared after the 
submission of our manuscript, we leave a detailed discussion of its 
results for future work. However, we note that in the particular 
simulations they presented, the cold classicals are dynamically ex- 
cited, inconsistent with the constraints we will establish. In Section 
15.1.31 we explore under what circumstances, if any, this could be 
avoided. 



sicals have a higher binary fraction because any hot clas- 
sical binaries were likely to have been disrupted when 
they were scattered from the inner disk to the classi- 
cal region (jParker fc Kavelaarsl 120101 ) and because bi- 
nary cap ture may have been less efficie nt in the in- 
ner disk (jMurray-Clay fc Schlichtingll201 If ) . However, to 
date this class of model has not yet been demonstrated 
to work quantitatively. We consider a generalized two- 
origin model in which hot classical deliver echoers as a 
result of scattering b y Neptune (r ather than due to reso- 
nance sweeping as in [Gomes 2003). Focusing on the con- 
sistency of this class of model with the eccentricity distri- 
bution of classical KBOs — unaccounted for by previous 
model realizations — we explore the parameter space for 
this generalized model using several alternative tactics: 

• Instead of attempting to produce a single model, 
we fully explore the parameter space of Neptune's 
eccentricity, semimajor axis, migration rate, eccen- 
tricity damping rate, and precession rate to assess 
the consistency of a collection of dynamical histo- 
ries with the observations. This approach is general 
in the sense that previou s models (e.g. Malhotra 
[19951: ILevison et all [20081 are under consideration 
(corresponding to a particular set of parameters), 
as well as other regions of parameter space that 
have not been explicitly considered. We will ex- 
plore Neptune's inclination and inclination damp- 
ing rate in a paper currently in preparation (R.I. 
Dawson and R. Murray-Clay 2012, in preparation). 
In Scction [5.4[ we clarify how to interpret complex 
solar system histories in the context of this general 
model. 

• Instead of matching the observations in detail, we 
focus on matching major qualitative features of the 
classical KBO eccentricity distribution that are un- 
affected by observational bias or by the long-term 
evolution of the solar system (i.e. the evolution 
that happens over the ~ 4 Gyr after the planets 
reach their final configuration). This approach al- 
lows us to perform short integrations that end once 
the planets reach their current configuration. 

• Instead of relying solely on numerical integrations, 
we determine which dynamical processes affect the 
evolution of the KBOs and place constraints using 
analytical expressions. 

• Instead of modeling all four planets directly, we 
model only Neptune but allow its orbit to change. 
We will demonstrate why this approach is sufficient 
for the problem we are exploring. 

Our exploration of parameter space could produce two 
possible outcomes. If we find regions of parameter space 
that can deliver the hot objects on top of the cold, these 
consistent regions will provide constraints for more de- 
tailed models. If we rule out all of parameter space, 
then a new type of model, employing different physi- 
cal processes, is necessary. Either way, we will identify 
and quantify what physical processes are responsible for 
sculpting the eccentricity distribution of the classicals in 
the generalized model we are treating. We emphasize 
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that, rather than proposing a new model, we are explor- 
ing a generalization of Neptune's dynamical history, in 
which previous models correspond to a particular set of 
parameters. 

In the next Section, we demonstrate that the hot and 
cold classicals have not only a bimodal inclination distri- 
bution, as already well established in the literature, but 
also distinct eccentricity distributions that were sculpted 
during Neptune's wild days. We use qualitative features 
of these eccentricity distributions to establish conserva- 
tive criteria that models must meet. In Section 3, we 
establish the framework for our study and argue that, 
combined with the distinct physical properties of the hot 
and cold populations, these eccentricity distributions im- 
ply separate origins for the cold and the hot classicals. 
In Section 4, we identify, for the classical region, the 
potential dynamical consequences of Neptune spending 
part of its dynamical history with high eccentricity - 
delivery of objects via scattering, secular forcing, acceler- 
ated secular forcing near resonances, and a chaotic sea — 
and present analytical expressions validated by numeri- 
cal integrations. In Section 5, we combine the analytical 
expressions from Section 4 with the conservative criteria 
established in Section 2 to place constraints on Neptune's 
path, orbital evolution timescales, and interactions with 
other planets, ruling out almost all of parameter space. 
We find that Neptune must spend time with high eccen- 
tricity to deliver the hot classicals by these processes, but 
is restricted to one of two regions of (a, e) space while its 
eccentricity is high. To avoid disrupting the cold objects, 
Neptune's eccentricity must have damped quickly or the 
planet's orbit must have precessed quickly while its ec- 
centricity was high. Finally, because Neptune's current 
semimajor axis is ruled out when Neptune's eccentricity 
is high, Neptune is constrained to have migrated a short 
distance after its eccentricity damped. In the final Sec- 
tion, we discuss our results and their implications for the 
early history of the solar system. 

2. CONSTRAINTS FROM THE OBSERVED ECCENTRICITY 
DISTRIBUTIONS OF HOT AND COLD CLASSICALS 

The cold and hot populations are defined by the ob- 
served bimodal inclination distribution of classical KB O s 
HBrowiJ2001t IGuibis et aLll20ToT: IVolk fc Malh"o1^l2lHll) . 
They also have distinct eccentricity distributions. The 
eccentricity distribution of all the observed KBOs is plot- 
ted in Figure [1] In Section 12. 1\ we present evidence for 
distinct hot classical and cold classical eccentricity dis- 
tributions and identify robust qualitative features of the 
distributions that models of Neptune's dynamical history 
must produce. In Section [2~2l we lay out the observational 
constraints which we will use for the remainder of the 
paper. In Section I2.3[ we assess the robustness of these 
features by performing statistical tests and considering 
observational bias. 

2.1. Evidence for distinct hot classical and cold classical 
eccentricity distributions 

We wish to use inclinations to separate the cold and 
hot classicals and then examine the eccentricity distri- 
butions of these two populations. Traditionally, the ob- 
served cold and hot objects have been separated using 
one inclination cutoff. However, because of the overlap 
between the hot and cold components in the bimodal 



CN fO fiT> 



" 0.2 




a(AU) 



Fig. 1. — Orbital eccentricity distribution of Kuiper Belt Ob- 
jects. The resonant and scattered objects are plotted as black 
pluses. The classical objects are plotted as colored circles. The 
red objects have i < 2° and are thus very likely cold classicals. 
The blue objects have i > 6° and are thus very likely hot clas- 
sicals. The membership of any givenpurple object, which has 
2° < i < 6°, is ambiguous (see Figure [2]l. Objects are taken from 
the Minor Planet Center Database and th e Canada-France Eclip - 
tic Plane Survey (CFEPS) and classified b y Gladman ct al. (2008); 
IKavelaars et alj !|2009fl ; IVolk fc Malhotral H201U 1. Dashed lines in- 
dicate the location of mean motion resonances with Neptune, which 
are included up through fourth order. 

inclination distribution, a single cutoff will necessarily 
result in the misclassification of hot objects as cold and 
vice versa. For example, if the classical population fol- 
lows the model KB O inclination distribution derived by 
IGulbis et al.1 ()2010f ) and we were to distinguish between 
the cold and hot populations using an inclination cut-off 
icut = 4°, 11% of objects with i < 4° would be actually 
be hot objects and 15% of objects with i > 4° would be 
actually be cold objects. Thus, using icut — 4° , 11% of 
the objects classified as cold would be "contaminated," 
and 15% of those classified as hot would be contami- 
nated. In Figure [H we plot the "contaminated" frac- 
tion over a range of values for icut for the cold and hot 
populations based on three models of the debiased in- 
clination distribution (iBrownl 120011: IGulbis et all 120101: 
IVolk fc Malhotral 1201 il ). For all three models, less than 
10% of the cold classicals are contaminated for icut < 2°, 
while less 3% of the hot classicals are contaminated for 
icut > 6°. 

Therefore, instead of using a single icut, we divide the 
classicals into a likely cold population (i < 2°), a likely 
hot population (i > 6°), and an ambiguous population 
(2° < i < 6°). We then examine the eccentricity dis- 
tributions of the likely cold and likely hot populations, 
which are "uncontaminated" samples. We use the un- 
contaminated eccentricity distributions to identify major 
features that models much match. In Appendix [5] we 
confirm that our results are consistent if we probabilisti- 
cally include the ambiguous population. 
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Fig. 2. — Fraction of each population "contaminated" by the 
other group as a function of the cut-off inclination icut between the 
cold and hot population. The dashed lines, labele d B01, are ca l- 
culated from the inclination distribution defined by Brown (2001); 
the dott ed lines, labeled G10 , from the inclination distribution de- 
fined by Gulbis ct al. (2010); and the dash-dot ted lines, labeled 
Vll, f rom the inclination distribution defined by Volk & Malhotra 

pour) . 

We wish to identify features of the eccentricity distri- 
bution that are sculpted during Neptune's wild days, not 
by the long-term stability of the region under the influ- 
ence of the modern solar system planetary configuration 
or by observational bias. First we compare the eccentrici- 
ties of observed likely cold (i < 2°) and likely hot (i > 6° ) 
objects to the survival map of Lykawka & Mukai (2005), 
.enera ted from a 4 Gyr simulation. Lykaw ka fe MukaH 
2005) generated initial conditions for test particles uni- 
formly filling a cube of (a, e, i) in the classical region: 
41.375.4t/ < a < 48.125 AU, < e < 0.3, and 
< i < 30°. They then performed a 4 Gyr numeri- 
cal integration including the test particles and the four 
giant planets (starting on their modern orbits). Then 
they computed the survival rate of test particles in bins 
of (a, e) and (a,i). In this work, we consider only the 
eccentricity survival map. This map bins over all incli- 
nationt0. After 4 Gyr of evolution under the influence of 
the planets in their current configuration, KBOs with an 
initially uniform eccentricity distribution would be dis- 
tributed according to this survival map. 

However, rather than following the survival map, the 
observed cold and hot objects exhibit major distinct 
features. In Figure we plot the sample of observed 
classical objects from the Minor Planet Center (MPC) 
and the Canada-Fran ce Ecliptic Plan e Survey (CFEPS; 
iKavelaars et al.1 120091 ) on top of the iLvkawka fe MukaH 
(|2005| ) stability map. The cold objects are confined to 
very low eccentricities. From 42.5 to 44 AU, the cold 
objects appear to be confined to e < 0.05. From 44 to 45 
AU, the cold objects appear confined below e < 0.1. This 
confinement of the cold classicals to below the survival 
limit implies that they were not excited above these levels 
because if they had been, we wo uld still observe objects 
at higher eccentricities. Similarly. Kavel aars et al.1 (2009) 
found that classical objects with i < 4.5° are restricted to 
42.5 AU < a < 45 AU and e < 0.1. In contrast, hot ob- 

3 We note that at a given semimajor axis, th e survival rate does 
not show a strong dependence on inclination iTCykawka fc Mukail 



2005, Figure 4, lower panel), except near the ug secular inclination 
resonance at 41.5 AU, which is devoid of low inclination objects 
We do not establish constraints in this region. 



jects occupy the upper portion of the survival region and 
appear uniformly distributed in a from 42 to 47.5 AU. 
Suggestively, they also appear to be distributed roughly 
along a scattering line, as if they were scattered into the 
classical region but did not have time to evolve to low 
eccentricities before Neptune's eccentricity damped. 

2.2. Conservative criteria that models must meet 

We use the following major qualitative features to place 
constraints on Neptune's dynamical history. We consider 
these criteria "conservative" because they allow for dy- 
namical histories at the very edge of consistency with the 
observations. 

Cold population: confined to low eccentricities of e < 
0.1 in the region from J^2.5 to 45 AU. In the region be- 
tween 42.5 and 45 AU, the cold objects have eccentrici- 
ties well below the distribution that follows the survival 
map. Therefore, Neptune cannot excite the cold classi- 
cal objects in this region above e = 0.1. (We choose this 
value to be conservative in ruling out reg i ons of param- 
eter space and to match [Kavelaars et al.l (j2009[ ) , but it 
appears that cold objects with semimajor axes less than 
44 AU are confined below e < 0.05, a tighter constraint.) 
We indicate this threshold as a solid yellow line in Figure 
1 

Hot population: delivered to the upper survival region 
with q > 34 AU out to 47.5 AU. The observed hot ob- 
jects occupy the upper portion of the survival region. 
Therefore, a consistent dynamical history should allow 
some objects to reach this region. It is not necessary for 
the transported objects to reach very low eccentricities, 
only low enough to survive under the current planetary 
configuration. We set the criterion that the hot classicals 
must be delivered to periapse q > 34 AU (dashed line in 
Figure [3]) from 42 A to 47.5 AU, the edge of observed 
population. 

2.3. Assessing the robustness of the observed features 

In determining which major features serve as con- 
straints on the dynamical history of the solar system, 
we address several complications: 

1. The inclinations of objects vary over time 
(IVolk fc Malhotrd[201Tl) . 

2. The inclination cut-off between the hot and cold 
classicals is model dependent. 

3. Proper elements are more robust than the observed 
instantaneous elements. 

4. The features in the eccentricity distributions might 
be the result of random chance or small number 
statistics. 

5. The eccentricity distributions may be impacted by 
observational bias. 

The first c ompli cation is addressed by 
iVolk fc Malhotral (|2011| ). They find that, at any 
given time, only 5% of objects will be inconsistent 
with their original inclination-based classification of hot 
versus cold. Therefore we expect the major qualitative 
features we identify to hold despite variations in the 
inclinations of some objects 
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Fig. 3. — Observed (plus symb ol) eccentricity dist r ibutio ns of cold classicals with i <° 2 (left) and hot classicals with i > 6° (right), 
plotted over the survival map of Lykawka & Mukai (2005). In the left panel, the cold classicals between 42.5-44 AU are confined to 
e < 0.05, well below the survival limit in this region, while cold classicals between 44-45 AU are confined to e < 0.1, also below the survival 
limit. In the right panel, the hot classicals occupy the upper portion of the survival region. We plot e = 0.1 from 42.5 to 45 AU as a 
solid yellow line in the left panel. The dashed li ne i n th e rig ht panel, periapse q = 34 AU, indicates an approximate upper threshold of 
long-term survival, which we will use in Scctions l5.2l and l5.3l Classical objects are taken f rom the Minor Planet Center Database a nd th e 
Canada-France Ecliptic Plane Survey (CFEPS; Kavclaars ct al. 2009) and are classified by Gladman ct al. (2008), Kavclaars ct al. (2009), 
and Volk & Malhotra (2011). 



The second complication is that the cut-off inclination 
between the hot and cold classicals depends on the pa- 
rameters and form of the model for the bimodal i nclina- 
tion distribution. The three models we conside r ([Brownl 
120011; IGulbis et aLll201d: iVolk fe Malhotrall201lD each use 
the functional form of sin i multiplied by a Gaussian but 
use different widths and cold/hot fractio ns. Th e y also 
use different planes for the inclination: I Brownl ([20011 ) 
defines the inclinat ion with respect to the ecliptic plane, 
IGulbis et al.l (|2010|) with respect to the m ean plane of the 
Kuiper Belt, and I Volk fc Malhotral (|201lD with respect to 
the invariable plane. However, despite these differences, 
i < 2° and i >° 6 are robust cut-offs for establishing an 
uncontaminated cold and hot population, respectively, 
for each of the three mo dels. In regards to the func- 
tional form of the model, IVolk fc Malhotral (|2011[) find 
that the high-inclin ation component is not we ll-described 
by a Gaussian, and Fabrvcky & Winn| (|2009f) argue that 
the most robust generic functional form for a distribu- 
tion of inclinations is a Fischer distribution. However, 
the discrepancies between different functional forms are 
strongest for classifying objects in the intermediate, over- 
lapping portion of the bimodal inclination distribution 
(K. Volk 2011, private communication), so we argue that 
our approach of definitively classifying only the "uncon- 
taminated" low and high inclination objects is robust. 

Regarding the th i rd com plication, the stability map of 
I Lykawka fe Mu kai (200]J) is formulated in terms of in- 
stantaneous eccentricity and inclination, but the most 



robust, no n- varying formulation of the orbital elements 
are the proper, or free, elements. However, none of the 
model inclination distributions are formulated in terms 
of the proper incl i nation , nor is the stability map of 
iLvkawka fc Mukai ([2005D formulated in the proper el- 
ements. To compare "apples to apples," we use the in- 
stantaneous orbital elements in the plots in this section. 
We use the proper elements in Appendix [Al and find that 
the observational features we identify (Section 12.21) still 
hold. 

Addressing the fourth complication, in Appendix [X[ 
we confirm that the confinement of the cold population 
to low eccentricities is statistically significant. For the 
hot population, we only impose the constraint that the 
models must deliver them to the long-term stable region 
(Section 12.21) ; we will demonstrate that this loosely for- 
mulated restriction ends up imposing strong constraints 
on Neptune's dynamical history. 

2.3.1. Ruling out observational bias through statistical tests 

We would not expect observational bias (the final com- 
plication) to cause the cold classicals to appear to con- 
fined to low eccentricities; KBO searches are more likely 
to preferentially observe high-eccentricity (i.e. small pe- 
riapse) objects of a given semimajor axis. However, to 
ensure that the features on which we base our constraints 
(Section 12. 2p are not created by observational bias, we 
perform the following test to see if observational bias 
could generate them: 
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1. We begin by generating a simulated sample of ob- 
jects uniformly distributed in (a, e). We set the in- 
clinations to follow the unbiased incli nation distri- 
bution of the classicals, as modeled bv lGulbis et al.l 

plch . 

2. Then we use the stability map of lLykawka fc Mukail 
(|2005f ) to transform this simulated sample fol- 
lowing a uniform eccentricity distribution into 
a sample following the eccentricity distribution 
shaped by the four giant planets under their cur- 
rent configuration. We call this process "filter- 
ing." For each simulated object, we obtain a pre- 
dicted survival rate from the (a, e) stability map of 
ILvkawka fc Mukail (|2005l ). Then we select a uni- 
form random number between and 1. If the 
randomly selected number is less than the pre- 
dicted survival rate, we include the object in our 
sample. Note that sin ce the survival rates of 
Lvkawka & Mukail (|2005|) are given as a range (e.g. 
10 - 20%, 90 - 100%), we repeated this entire 
test (i.e., steps 2-4) three times, once using the 
minimum of each range, once using the mean of 
each range, and once using the maximum of each 
range. As expected, the resulting eccentricity dis- 
tribution had higher (lower) eccentricities when we 
used the maximum (minimum) each range but the 
major features we identified still held. The simu- 
lated population in Figure 0] uses the mean. Next 
we transformed the "survival-rate filtered" sample 
from step 2 into an observed sample: 

3. We randomly assign each object an H magnitude^ 
between 6 and 8. 

4. Then we apply the L7 Su rvey Simulator for th e 
well-characterized CFEPS. (|Kavelaars et al.ll2009D . 

5. We compare the final simulated distribution to the 
subset of objects tha t were detected by CFEPS 
(jKavelaars et al.ll2009l ) (FigureE}. Figure 0] is anal- 
ogous to Figure El It includes the simulated distri- 
bution (circles) and only the subset of KBOs ob- 
served by CFEPS. The simulated distribution is 
not confined to e < 0.1 from 42.5 to 45 AU, con- 
firming this feature of the observed eccentricity dis- 
tribution does not result from observational bias. 
Note also that simulated hot objects are found at 
lower eccentricities than observed. 

3. FRAMEWORK 

Our representation of the observations in Figure |3] 
highlights the problem with theories of a single origin for 
the hot and cold objects. How could a single origin pro- 
duce both a cold population confined to low eccentricities 
and a hot population, with different physical properties 
and inclinations, dwelling above at high eccentricities? 
In Section 13.11 we explain why a single origin scenario is 
unlikely In Section [3.21 we describe a scenario, which we 
will consider throughout the rest of the paper, in which 

4 An alternative method, randomly drawing the H magnitudes 
from observed classical CFEPS objects, yielded results that were 
qualitatively the same. 



the cold population formed in situ and the hot popu- 
lation formed in the inner disk and was transported to 
the classical region. In Section 13.31 we explain why is 
it reasonable to place constraints on Neptune's history 
using the evolution of the KBOs only during Neptune's 
wild days, and in Section 13.41 we discuss the possibility 
of alternative scenarios of Kuiper Belt assembly. 

3.1. Ruling out a single origin for the hot and cold 
classicals 

A single origin for the hot and cold populations seems 
unlikely If the cold and hot classicals formed together 
in the classical region, where they are observed today, 
it is difficult to imagine a process that would excite the 
hot population while leaving the cold populatio n confined 
to low eccentricities. lHahn fc Malhotral (|2005l ) proposed 
a scenario in which the classical region has been pre- 
excited. However, this scenario does not produce a pop- 
ulation of cold classicals confined to low eccentricities. 
Moreover, if both the cold and hot classicals were trans- 
ported from the inner disk, it seems unlikely that a com- 
mon deposition process wo uld place the cold classicals 
sole ly at low eccentriciti es. iLevison fe Morbidellil (2003) 
an d lLevison et~a l. (2008) propose scenarios in which both 
the hot and cold classicals are transported from the inner 
disk. 

In the scenario of ILevison fe Morbidellil J2003), the 
cold classicals were pushed outward by the 2:1 resonance 
and dropped during stochastic migration, while the hot 
classicals scattered off of Neptune. The feasibility of this 
mechanism depends on the size distribution of planetesi- 
mals, because the migration needs to be stochastic in or- 
der to drop objects from resonance. When Neptune scat- 
ters a planetesimal inward and the planetesimal is ejected 
by Jupiter, Neptune experiences a net gain in angular 
momentum and migrates outward. If the planetesimals 
are small, this is a smooth process, but if they are large, 
it is a jumpy, stochastic proce ss, in which KBOs can 
be dro pped from resonance. See iMurrav-Clav fc Chiang! 
(2006) for a detailed exploration of stochastic migration; 
they conclude that planetesimal-driven migration cannot 
generate the necessary stochasticity unless a large frac- 
tion of planetesimals formed very large. This constraint 
merits a fresh lo ok in light of new plan etesimal forma- 
tion models (see lChiang fc Youdinl l2010l and references 
therein). However, even if extreme planetesimal prop- 
erties allowed this mechanism to work, objects dropped 
from the 2:1 resonance would have a range eccentricities, 
not be confined solely to low eccentricities. Therefore 
this mechanism holds more potential for producing the 
hot population than the cold popu l ation . 

In the scenario of ILevison et al . (2008), the cold clas- 
sicals are objects that, like the hot classicals, were scat- 
tered into the Kuiper Belt by an eccentric Neptune but, 
unlike the hot classicals, evolve down to low eccentrici- 
ties in regions near resonances. However, this mechanism 
would create a range of eccentricities for the cold classi- 
cals and thus would have trouble producing the confined 
eccentricities of the region of 42.5 AU < a < 45 AU 
(Figure [3]). They do find some correlation between a 
particle's final inclination a nd its initial semima jor axis 
in one of their simulations (ILevison et al.l 120081 Figure 
11, panel (b)), which may be able to partially account 
for a difference in physical properties between low and 
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Fig. 4. — CFEPS objects only. Plotted over the survival map of Lvkawka & Mukai (2005) are predicted (circle) and observed (plus) 
distributions of cold classicals with i < 2° (left) and hot classicals with i > 6° (right). The predicted classicals are t he distribution expected 
from a uniform (a, e) distribution, filtered by the survival map and put through the CFEPS Survey Simulator of Kavclaars ct al. (2009) 
(see the text for further detail). The observed classicals are those observed by CFEPS. In the left panel, the cold classicals between 42.5 
and 44 AU are confined to e < 0.05, well below the survival limit, while cold classicals between 44 and 45 AU are confined to e < 0.1, also 
below the survival limit in this region. We plot e = 0.1 from 42.5 to 45 AU as a solid yellow line in the left panel. In the right panel, 
the hot classicals occupy the upper po rtion of the survival region. The dashed line indicates an approximate upper threshold of long-term 
survival, which we will use in Sections 15.21 and 15.31 



high inclination objects. However, transporting the cold 
classicals from the inner disk is not consistent with the 
finding by Parker & Kavclaars (2010) that wide binaries 
- of which the cold population contains a number - 
cannot survive transportation from the inner disk to the 
classical region. 

3.1.1. The eccentricity distribution was not sculpted solely 
by a different stability threshold in the past 

One might wonder whether the observed confinement 
of the cold classicals (Figure [3]) is the result of a smaller 
stability region than exists today, as if the cold classicals 
follow an ancient scattering line. However, there are nu- 
merous "hot" objects with i > 6°, as well as ambiguous 
objects with 2° < i < 6°, in the region from 42.5 to 45 
AU that have high eccentricities, right up to the modern 
stability limit. To create the observed distribution, one 
would need a mechanism that removes all objects with 
high eccentricities and i < 2° while leaving a) objects 
with low eccentricities and i < 2°, and b) objects with 
a range of eccentricities and i > 2°. Therefore, it seems 
unlikely that this mechanism could produce both the hot 
and cold populations. We note that such a scenario could 
take place before Neptune transports the hot classicals. 
However, such initial sculpting would not affect the con- 
straints we will place, which Neptune still needs to obey 
during the hot classical transport phase. 

3.2. Colds in situ, hots transported from the inner disk 



Thus, throughout the rest of paper, we consider the 
genera l scenario — also discussed in iMorbidelli et al.1 
(2008) — in which the hot objects are transported to 
the classical region from the inner disk and the cold ob- 
jects form in situ in the classical region. The cold objects 
must not be dynamically excited, as quantified by the 
criterion we established in Section 12.21 In Figure we 
show a conceptualization of this model. This general sce- 
nario encompasses previous models and allows Neptune 
to undergo any potential combination of high eccentric- 
ity, migration, and/or apsidal precession with a range 
of initial eccentricities and semimajor axes. If our con- 
straints do not rule out all of parameter space for this 
model, it may be possible to produce both the hot and 
cold classical population. Otherwise, a major physical 
process is missing from current models. 

For now, we can think of Neptune as having a high ec- 
centricity at one location after undergoing planet-planet 
scattering. Our results hold in more complicated scenar- 
ios as well, as we will describe in Section 15.41 



3.3. The case for considering short-term evolution 

We have chosen features (Section 12. 2|) that are not 
shaped by the long-term survival of KBOs under the cur- 
rent solar system planetary configuration. Therefore, we 
can focus on modeling the processes that affect these fea- 
tures during the interval of Neptune's wild days instead 
of treating the entire 4 Gyr. We model these processes 
analytically in Section [4] and validate our analytical ex- 
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Fig. 5. — Conceptual framework, (i) The hot classicals (blue) form in the inner disk and the cold classicals (red) form in the outer disk, 
(ii) Neptune is scattered onto a high-eccentricity orbit through its interactions with other planets, (iii) Neptune scatters the hot objects 
into the classical region without disrupting the cold ones, (iv) Neptune's large eccentricity damps, leaving the cold classicals confined to 
low eccentricities. 

pressions using numerical integrations. 

Unless otherwise specified, the integrations are per- 
formed as follows. We perform the integra tions using 
the M ercury 6.2 hybrid symplectic integrator (Chambers 
1999), an iV-body code that allows massless test parti- 
cles. We employ an accuracy parameter of 10 -12 and 
a step size of 200 days and impose (if applicable) the 
migration and damping of Neptune's orbit through user- 
defined forces and velocities. Migration and damping 
follow the following forms: 

a N = (o N )/ + ((on)o - (oN)/)exp(-t/r ON ) 
eN = (e N ) exp(-t/T eN ), 

(1) 

where <2n is the semimajor axis of Neptune at time t, 
(on)o is the initial semimajor axis, and (on)/ = 30.1 
AU is the final semimajor axis. At time t, the eccen- 
tricity of Neptune is e-^; Neptune's initial eccentricity is 
(eN)o- The forced evolution of Neptune's orbit is imple- 
mented through modifications to Mercury 6 .2, de scribed 
in detail in the Appendix of iWolff et all (|2012t ). The 
migration and damping are parametrized by timescales 
r aN and r eN respectively, which we specify in the text in 
the applicable cases. When specified, Neptune is forced 
to undergo apsidal precession using an artificial stellar 
oblateness force, built into Mercury 6.2, that we mod- 
ified to apply only to Neptune, parameterized by a Ji 
coefficient chosen to produce the correct precession rate. 
Unless otherwise noted, Neptune is the only planet in- 
cluded in the integration. The KBOs are modeled as 
600 massless test particles with initial a evenly spaced 
between 40 and 60 AU and initial e = i = 0. The migra- 
tion, damping, and apsidal precession are not applied to 
the KBOs, only to Neptune. 

3.4. Alternative scenarios for Kuiper Belt assembly 

The purpose of this paper is to explore whether the 
generalized scenario described here can even work, i.e. 
whether it is ever possible to transport the hot classi- 
cals from the inner disk to the classical region without 
disrupting an in situ cold population. Yet alternative 
scenarios exist that do not fit within this framework, 
such as the additional pla net beyond Pluto proposed by 
Lykawka & Mukai (2008) and others. Furthermore, we 
will describe in the conclusion how additional constraints 
could rule out the generalized model we consider. In that 
case, development of alternative scenarios would be nec- 
essary. Obviously, the constraints we will place do not 
necessary hold for a scenario that is not encompassed by 
our general model. 

4. PHYSICAL PROCESSES RESULTING FROM NEPTUNE'S 
HIGH ECCENTRICITY 



We begin our analysis by describing the physical pro- 
cesses that can impact the Kuiper Belt if Neptune's ec- 
centricity is high. First we consider how classical KBOs 
reach the classical region. In our generalized model (Fig- 
ure [5J , the cold classicals form in situ and hot classi- 
cals are delivered by Neptune from the inner disk. Once 
the hot and cold objects are in the classical region, both 
evolve in response to an eccentric Neptune. A tension 
arises between need of hot objects for Neptune to be ec- 
centric — to deliver them into the classical region and to 
allow them to evolve to low eccentricities once they arrive 
- and the undesirable excitation of cold objects when 
Neptune is eccentric. In this Section, we lay out analyt- 
ical expressions for how KBOs evolve in response to an 
eccentric Neptune and use this theory to transform our 
observational constraints into comprehensive constraints 
on Neptune's orbit during its high eccentricity period. 
We will employ these constraints on Neptune's orbit in 
Section [5] to rule out much of parameter space. 

4.1. Delivery into the classical region 

In the generalized model we explore, Neptune may be 
scattered outward from the inner solar system onto a 
highly eccentric orbit. After this occurs, Neptune's new 
orbit crosses the orbits of some planetesimals in the inner 
disk (see Figure [5]), which scatter off the planet. This 
mechanism can potentially deliver hot objects from the 
inner disk into the classical region. 

The region into which Neptune can scatter objects is 
defined by the planet's semimajor axis on and eccentric- 
ity eN- Neptune can scatter objects outward to peri- 
apses q — a(l — e) between Neptune's periapse r p ^ = 
(on - a H )(l — eN) and apoapse r a ^ = (a N +a H )(l + e^). 
In Figure El we show examples of the region into which 
Neptune can scatter KBOs for two sets of parameters 
(<ZN,eN). We have adjusted r Pi ^ and r Qj N to include 
the Hill sphere radius, an (~ 1 AU), the distance from 
Neptune at which Neptune's gravity overcomes the Sun's 
tidal gravity. Particles that enter Neptune's Hill sphere 
will be scattered and, if they are scattered outward into 
the classical region, will reach - under the approximation 
that the scattering location becomes the particle's new 
periapse - a given semimajor axis a with eccentricities 
between 1 — r p ^/a and 1 — r a ^/a. 

4.1.1. Implications of scattering for the cold classicals 

If Neptune's apoapse r a ^ is large enough, the planet 
can potentially impinge into the cold classical region, 
scattering the objects there. However, we find in prac- 
tice that if Neptune's eccentricity damps within the con- 
straints we will place, then cold objects are never ex- 
cited up to the scattering line, which moves upward as 
Neptune's eccentricity damps. In Section 15.1. 31 we will 
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Fig. 6. — Region into which Neptune can scatter particles for 
apj = 24, eN = 0.02 (top), a typical initial condition for extensive 
migration models, and ajj = 2 8.9, = 0-3 (bottom ), which is the 
initial condition for Run B in[Levison ct al. (20QH). The dashed 
line marks the upper threshold of long-term survival, as indicated 
in Figure [4] 

return to this issue for the special case of Neptune un- 
dergoing fast apsidal precession. 

4.2. Secular forcing 

Secular forcing has a strong effect on the behavior of 
KBOs during the poten tial period when Neptune 's ec- 
centricity is high. (See iMurrav fc DermottJ 120001 for a 
pedagogical description of first order secular theory out- 
side of mean motion resonance.) On a timescale of a 
million years, a cold object is excited to a higher eccen- 
tricity via Neptune's secular forcing. (The direct forcing 
from the other planets on the KBO is negligible and the 
other planets only affect the KBOs via Neptune, as we 
will demonstrate in Section f4. 2. 3p . A hot object that is 
scattered into the classical region also experiences sec- 
ular forcing, which can decrease its eccentricity so that 
its orbit no longer crosses Neptune's. Hot or cold, an 
object's eccentricity is a vector combination of its forced 
eccentricity — set by Neptune's eccentricity, Neptune's 
scmimajor axis relative to the object's, and Neptune's 
apsidal precession rate — and of the object's free eccen- 
tricity, which is set by its initial condition before Nep- 
tune is scattered to a high eccentricity. The object's 
free eccentricity precesses about the forced eccentricity 
at the secular frequency <?kbo- Therefore, as we will 
demonstrate, a cold object has a well-defined excitation 
time and amplitude, and a hot object will have a mini- 
mum eccentricity it can reach after being scattered into 
the classical region. Thus while Neptune's eccentricity is 
high, secular forcing potentially is an important mech- 
anism for exciting the cold objects and stabilizing the 
hot ones. As we will show, when Neptune's eccentricity 
damps quickly, the orbits of the KBOs are "frozen" near 
the eccentricities they reached through secular evolution. 

4.2.1. Basic secular evolution 

First we define expressions for the secular evolution of 
a test particle under the influence of a planetary system 
containing only Neptune and the Sun. The components 
of the particle's eccentricity vector are h = e sin w and 
k = ecosTu, where tu is the particle's longitude of pe- 
riapse. Secular forcing by Neptune causes h and k to 
evolve as (to first order in e and e^): 

h = ef re e sin^KBO* + P) + eforced sui(tun) 



(2) 
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The constants ef ree and /3 are determined from the ini- 
tial conditions, and the particle's forced eccentricity is 
etorccd- Here, tun is the longitude of periapse of Nep- 
tune, eN is the eccentricity of Neptune, and a is the 
ratio of Neptune's semimajor axis to that of the particle, 
all of which are assumed to be const ant. The functions b 
are st andard Laplace coefficients (see lMurray fc DermottJ 
2000). The secular frequency of the KBO is <?kbo, wn 
is the mass of Neptune, is the mass of the Sun, 

n = (Gm^/a 3 ) 1 / 2 is the particle's mean motion, and G 

is the universal gravitational constant. 

Consider a cold object with e = at t = 0. In our 
approximation, Neptune, having been scattered by the 
other giant planets, effectively instantaneously appears 
and imparts a forced eccentricity of ef orcc d- From these 
initial conditions, ef rcc = ef orcc d. Then, the KBO's forced 
eccentricity vector remains fixed and the object's total 
free eccentricity vector precesses about the forced eccen- 
tricity. Thus its total eccentricity varies sinusoidally from 
e = to e = 2ef orco d on a timescale set by <7kbo- 

A hot object scattered into the classical region, in con- 
trast, has an eccentricity e at t = 0. The magnitude of 
its free eccentricity is a value between max(0, e — ef orce d) 
and e + ef orce d, depending on the initial location of its 
periapse relative to Neptune's. Over a timescale set by 
<?KBO; its total eccentricity oscillates. Depending on the 
initial conditions, it may reach an eccentricity low enough 
so that its orbit no longer crosses Neptune's and/or so 
that it is stable under the current configuration of the 
giant planets. 

An example of the secular evolution of cold objects 
"going up" and hot objects "going down" in eccentricity 
is shown in Figure [7J highlighting the tension between 
the evolution of hot objects to low eccentricities and the 
evolution of cold objects to high eccentricities. The cold 
objects (red) begin with e = (see Section 13.31 for a 
general description of the integrations we performed.) 
The hot objects in the integration (blue) all begin with 
e = 0.2 and tu = tun + 7r/3, for the purposes of illus- 
trating secular evolution^- As time progresses through 
three snapshots, the cold objects become excited and 
the hot objects reach low eccentricities. The analytical 
model from Equation @ matches well except near mean 
motion resonances, where the secular evolution is much 

5 As shown in Scction l4.ll a real hot object can only be scattered 
to a certain region of (a, e) space in the classical region, and its 
eccentricity and periapse are actually correlated 
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faster than predicted. We also overplot a more accurate 
analytical expression that includes a resonant correction 
term, which we will derive in Section [4.2.21 

We note that throughout the paper, we perform numer- 
ical integrations of objects with initial semimajor axes 
out to 60 AU to give a better conceptual picture of the 
secular excitation. Moreover, depending on where the 
initial population was truncated in the solar system's 
planctesimal disk, it is possible that additional classi- 
cal KBOs will be discovered beyond 48 AU in the future, 
and we would like to make testable predictions. Finally, 
the integration results should be interpreted as examples: 
since we cannot show a figure for every possible combina- 
tion of parameters for Neptune, we plot out to 60 AU to 
let the reader imagine the results if Neptune's semimajor 
axis were smaller. 

4.2.2. Refined secular expression 

The secular expression in Section l4.2.1l is valid to first 
order in eccentricity, neglects the effects of orbital reso- 
nances, and applies in the case in which Neptune's orbit 
does not apsidally precess. However, these neglected ef- 
fects can significantly alter a secularly evolving KBO's 
behavior: 

1. Neptune's high eccentricity makes terms of order 
ejjj (ignored in deriving Equation^ non-negligible. 
As we will show, these extra terms result in a faster 
secular forcing frequency of the KBO. 

2. Proximity to mean- motion resonances with 
Neptune significantly alters the secular fre- 
quencies of KBOs (as s hown f or other solar 
syste m bodies in iHilll Il897t iMalhotra et al.l 
1989: iMinton fc Malhotral 120111). Following 
Malhotra et al.l (|1989) , we incorporate resonance 
correc tion terms, described in our Appendix 
IB. 1.21 These resonance correction terms are very 
important for objects near resonance but not 
valid for objects librating in resonance (we are not 
considering resonant objects in this papei0). 

3. As Neptune's eccentricity is damped, the particle's 
forced eccentricity goes to zero. If the damping oc- 
curs over a timescale r eN shorter than the secular 
oscillation time, the particle's eccentricity is frozen 
at approximately the value it reaches at the eccen- 
tricity damping time. If the damping occurs over a 
longer timescale, the particle's eccentricity evolves 
to its initial free eccentricity. 

4. Apsidal precession of Neptune alters the forced ec- 
centricity, keeping it low when Neptune precesses 
quickly. It also alters the oscillation timescale of 
the total eccentricity, because now both the free 
and forced eccentricity are precessing. 

5. Migration alters the secular frequencies and shifts 
the locations of the resonances. 

6 It has been claimed (e.g. Lcvison et al. 2008) that the entire 
classical region from the 3:2 to the 2:1 re sona nce is full over over- 
lapping resonances. However, in Section 14.51 we will demonstrate 
that Neptune's high eccentricity alone does not cause resonance 
overlap in the classical region. 



Based on these considerations, we modify the standard 
expression (Equation [2]) for the secular evolution of a 
test particle's h and k under the influence of an eccentric 
Neptune to incorporate these effects. 

h = ef ree ,0 Sm(g KB ot + A)) + ^forced sin(a7 N , + tn N t) 
k = e frcC! o cos(g K Bo£ + A)) + forced cos(tu N . + ro N t) 

(4) 



5kbo = (1 + —e N )ab\/ (a) + ed.g K BO (5) 

^KBO = a(Cxe^ 1 ) 2 ^n(6) 
m 

O 

eforcod = sin(min(g gKBO T eN ,7r/2)) , Skbo e N (t) (7) 

— <?KBO 

9kbo = - 1 + 7~e N ) ab\ L(a) (8 

fix °> m 4 

O 

Compare Equation (j4]) to Equation ([2]). The form is 
the same but Equation Q has several key differences 
and new variables, that we will proceed to discuss and 
define throughout the remainder of this subsection. One 
distinction is that several quantities that were fixed in 
Equation ^fy (i.e., w^,a, ejsr) can now vary with time. 
Additionally, we have incorporated corrections for Nep- 
tune's high eccentricity and for the potential proximity of 
the KBO to mean motion resonance with Neptune. The 
impact of the time- varying eN is the most complicated, 
so we leave it for last. 

As our first correction, we allow the longitude of peri- 
center of Neptune, zjjn, to precess. Neptune's precession 
adds an extra term w-^ to Equation ^ (where we rewrite 
ton as a linear function of time: wn.o + ^n4) and to 
Equation (J7J, in analogy to the standard secular theory 
for the four-planet case. We note that the other giant 
planets impact the secular evolution of the KBOs indi- 
rectly by causing Neptune's eccentricity to precess (see 
Section 14.2.31 for discussion) . 

The precession rate sjn has the same role (and same 
place, in the denominator of the forced eccentricity) in 
the just-Neptune secular theory as in t he standard four- 
planet secular theory (Appendix IB.2|) for the particle 
except that we are specifying Neptune's evolution via 
w'n instead of constructing a secular theory for the plan- 
ets that produces a particular cJn- When is large, 
the effective forced eccentricity ef orco d = |ef or cod|, de- 
fined below, remains low because the forced eccentric- 
ity is inversely proportional to the precession rate for 
I^n| > I^kboI (Equation 13) . 

When Ne p tune migrates, a changes with time. In 
IWolff et all (l20TD . we found that migration occurs 
in three regimes, relative to the eccentricity damping 
timescale r eN : fast, comparable, and slow. When Nep- 
tune's migration timescale r aN is slow relative to the 
damping time, the secular evolution of the KBOs effec- 
tively takes place as if Neptune remains at its initial lo- 
cation. When Neptune's migration is fast relative to the 
damping time, the secular evolution of the KBOs effec- 
tively takes place as if Neptune were always at its final 
location. In the intermediate case, in which r eN ~ t Qn , 
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Fig. 7. — Numerical integration (an = 30, = 0.2) shows cold objects (red) secularly evolving to high eccentricities and hot objects 
(blue) to low eccentricities. The gray is our analytical model without the resonance correction terms (Equation |2j and the black, which 
matches much better near the resonances, includes resonance terms (Equation [4} . 



modeling the secular evolution at the location Neptune 
reaches after half a damping time is a fair approximation. 
In this work, we therefore use a fixed a, which should be 
chosen according to these principles for a given evolution 
of Neptune. Figure which we will describe after dis- 
cussing our treatment of Neptune's eccentricity, provides 
an example showing that this approach is effective. 

Before considering the impact of eN varying with time, 
we discuss the correction terms for Neptune's high ec- 
centricity and for resonances. The correction terms for 
Neptune's high eccentricity do not change the form of 
the secular evolution. We have applied a correction term, 
■^■ejSj, in our expression for <?kbo (Equation [5]) , and an- 
other such factor, ^-e^, in g' KBO (Equation [8]) , which 
is another eigenfrequency. These terms are derived, the 
/ factors (which are of order unity) defined, and their 
necessity demonstrated in Appendix [Bj 

Proximity to resonance changes the se cular frequency 
qkbo (Equation [SJ, as described by iMalhotra et al.l 
(1989). Orbital resonances greatly increase the secular 
forcing frequency because terms in the disturbing func- 
tion that depend on the resonant angle can no longer be 
averaged over. The amplitude of the resonant correction 
term e is defined in Appendix |B] and depends on how 
close the particle is to the location of mean-motion res- 
onance. The frequency, 5(7kbo, defined in Equation ©, 
depends on the order of the resonance x and a constant 
C, of order unity, that is different for each resonance. See 
Figure [7] for a demonstration of the resonance correction 
terms. 

Finally, we turn to the impact of eccentricity damp- 
ing, which alters eN- Instantaneously, the KBO 



has the eccentricity components h = e sin w = 
ef r ee sin <f> + eforced sin ton and k — e cos zu — ef ree cos 4> + 
eforcod cos wn, where 4> = <7kbo t + j3. However, now the 
forced eccentricity vector (ef orce dCOSCi7N, ef orce d sintUN) is 
changing. We have already accounted for the apsidal 
precession, but now the magnitude of the forced eccen- 
tricity vector is changing as well. When ef orcc d does not 
change, ef rcc is a constant determined by initial condi- 
tions. This remains true when ef or ccd evolves slowly com- 
pared to the secular forcing time of the KBO. Otherwise, 
ef rco changes. Instead of allowing both ef ree and ef or cod to 
change with time, we use ef rec ,o and define an "effective" 
forced eccentricity, |ef orce d| (Equation[7]). The "effective" 
forced eccentricity changes with time proportionally to 
cn- Throughout the rest of the paper, we will refer to 
efrcc. o as efroo and the "effective" forced eccentricity as 
eforced- The "effective" forced eccentricity includes a fac- 
tor sin(min((? ffKBO T eN , tt/2)). When Neptune's eccentric- 
ity damping timescale is long compared to the KBO's 
secular oscillation period (i.e. jkbo^n > 7r/2) , the par- 
ticle's total eccentricity damps to its initial free eccen- 
tricity, ef rcc ,o, and this factor is unity. However, when 
Neptune's eccentricity damps quickly, the particle's to- 
tal eccentricity damps to a value near the eccentricity 
it reached after one damping time. The empirical cor- 
rection factor allows us to model the particle's behavior 
without altering the form of the secular evolution. This 
empirical factor provides a match to the integrations (see 
Figure E|). 

The modified secular theory, Equation Q, matches 
the integrations even when damping and migration are 
included (Figure [8]). Figure [8] shows an example of a 
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case in which r eN = 0.3 Myr is shorter than its migration 
timescale, r aN = 5 Myr. (See Section [3731 for a descrip- 
tion of how we implemented the damping and migration 
of Neptune's orbit.) In the top row, Neptune undergoes 
eccentricity damping from e = 0.3 at constant semima- 
jor axis 28 AU. In the middle row, Neptune migrates 
from 28 AU to 30 AU on a timescale t Qn = 5 Myr, with- 
out its eccentricity damping. Unlike any other model in 
this paper, we use a time-dependent a(t) for the analyt- 
ical model. In the bottom row, Neptune undergoes both 
damping and migration, but we plot again the same an- 
alytical model as in row 1, in which Neptune undergoes 
only eccentricity damping (no migration). The model 
over plotted in row 3, even though it in no way includes 
the effects of migration, matches well. When the migra- 
tion is slow compared to the damping, the change in the 
secular frequency Qkbo is negligible over the timescale 
during which Neptune's eccentricity is high. It is as if 
Neptune's eccentricity damps while Neptune remains at 
its initial on- Thus when the migration is slow compared 
to the damping timescale, we can model the KBOs' sec- 
ular evolution as if Neptune's eccentricity damps while 
Neptune remains in place. 

4.2.3. Effects of other planets 

Our approach of modeling only Neptune, undergoing a 
range of orbital histories that would be caused by in- 
teractions with the other giant planets and with the 
solar system's planetesimal disk, is sufficient because 
the other planets primarily affect the KBOs only indi- 
rectly, through influencing Neptune. The influence on 
the Kuiper Belt of a single, apsidally precessing Neptune 
matches the influence of multiple planets in both inte- 
grations and theory. 

An illustrative case is shown in Figure The initial 
conditions for the particles are the same as in Figure 8, 
and Neptune has on = 28 and eN = 0.2. The top panel 
(four planets) and middle panel (just Neptune undergo- 
ing precession) are very similar, keeping the objects at 
lower eccentricities than in the bottom panel (just Nep- 
tune, no precession). Thus precession must be included 
in the parameter space exploration, and including preces- 
sion successfully accounts for the influence of the other 
giant planets. We note that in the time of the snapshot 
(1.4 Myr), the particles have secularly evolved to have ec- 
centricities large enough so that their orbit, at the proper 
orientation, could intersect Neptune's. In cases in which 
Neptune precesses (top two panels), the objects are scat- 
tered by Neptune as Neptune's orbit precesses to inter- 
sect the orbits of the KBOs with eccentricities above the 
scattering line. The cut-off is at 45 AU because, interior 
this location, particles are secularly evolving quickly due 
to their proximity to resonance and have thus reached 
high eccentricities, allowing them to scatter. In the final 
panel, the particles are not scattered because Neptune's 
orbit does not precess to intersect the orbits of the par- 
ticles. 

The other planets matter in that they affect Neptune, 
but their direct effect on the KBOs is negligible. Their 
main effect is to cause precession of Neptune. Interac- 
tions between Uranus and Neptune, if they closely ap- 
proach the 2:1 resonance while Neptune's eccentricity is 
high, cause additional effects that we discuss in Section 
14.51 These effects can also be modeled using only Nep- 



tune, with appropriate orbital variations. 

Treating only Neptune reduces the number of free pa- 
rameters, allowing a more thorough exploration of the 
restricted space. The constraints we will develop can be 
applied to more extensive m odels that include the other 
giant planets. See Appendix IB. 21 for a mathematical dis- 
cussion of how the full four-planet secular theory reduces 
to the just-Neptune case. 

4.2.4. Constraints from secular excitation of cold objects 

It follows from the expressions in Section B. 2.21 that the 
excitation of the cold classicals happens on timescales 
of millions of years (or shorter near resonances) with an 
amplitude and timescale that depend on Neptune's semi- 
major axis, eccentricity, eccentricity damping timescale, 
migration rate, and precession rate. Complementarily, 
hot objects scattered into the classical region can evolve 
to lower eccentricities on similar timescales. 

The observations require that the cold classicals not be 
excited above e > 0.1 in the region 42.5 AU < a < 45 
AU, as demonstrated in Section [21 From Equations Q 
and ([5]) derived in Section 14.21 it follows that the cold 
classicals will not be excited above e > 0.1 at a given 
location if (CONSTRAINT 1): 



sin(min(g K BOT eN ,7r/2))l , |e N < 0.1 (9) 

»N — ffKBO 

Thus there are three possible regimes to "preserve" an 
in situ population of cold classicals through Neptune's 
wild days: 

1. The eccentricity of Neptune is small enough that 
the region's forced eccentricity, proportional to eN , 

is below 0.1 (i.e. I . 9 ' KBO |e N < 0.1) 

2. Neptune's periapse precesses quickly enough that 
the region's forced eccentricity, inversely propor- 
tional to ton, is below 0.1 (i.e. ton is large). 

3. The eccentricity of Neptune damps quickly enough 
that the objects are not excited above 0.1 (i.e. T e 
is small). 

4.3. Effects of post-scattering secular evolution on hot 

objects 

Hot objects that have been scattered into the classical 
region from the inner disk (Section l4.1[) will undergo sec- 
ular evolution when they arrive in the classical region. 
They will reach a given semimajor axis a in the clas- 
sical region with eccentricities between 1 — r a ^/a and 
1 — r p n/ci. Not all of these eccentricities are consistent 
with stable orbits over 4 Gyr (Figure [3]) . If a particle 
is scattered to a high eccentricity above the stable re- 
gion, under certain conditions — if T eN is not too fast 
and Neptune imparts a forced eccentricity that is large 
enough relative to the particle's free eccentricity — the 
particle can reach a region of long-term survival through 
secular evolution. 

In Figure [TU1 we show two examples of KBOs that are 
scattere d into the classical reg ion from an integration re- 
sembling lLevison et al.l (|2008|) Run B. In this integration, 
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Fig. 8. — Secular excitation of cold classicals when Neptune undergoes eccentricity damping and/or migration. In the top row, Neptune's 
eccentricity damps with r eN = 0.3 Myr, and the planet does not migrate; secular theory, including damping (Equation [4{ , is plotted in 
black. In the middle row, Neptune migrates outward on the timescale Ta N = 5 Myr, and its eccentricity does not damp. The secular theory, 
including migration (Equation [4j, is plotted in black. The scattered points in the final panel are objects that have undergone accelerated 
forcing near resonances as the resonances swept through. The bottom row displays a numerical integration including both eccentricity 
damping and migration, on the same timescales as above, but the analytical model overplotted in black is the same as in row 1. 

panel) but secularly evolve down into the stable region. 
After 4 Gyr, the particles remain in this location. Nei- 
ther of these example objects is librating in an orbital 
resonance. 

Therefore, Neptune's apoapse r p ^ must be large 
enough so that the hot classicals reach the region of long- 
term survival either immediately or can evolve there be- 
fore the eccentricity of Neptune damps, after which the 
eccentricity of the particle is frozen. 

However, the post-scattering evolution depends 
strongly on vo — w^- Particles are not actually scattered 
to orbits with independent, random w — w^. This is be- 
cause, by definition, after each scattering, the particle's 
new orbit fulfills the condition that at 9, the angle at 
which the orbit of the particle and the orbit of Neptune 



Neptune begins with ajv = 28.9 AU, eN = 0.3. Its eccen- 
tricity damps on a timescale of r e = 2 Myr, and it un- 
dergoes migration to 30.1 AU with r Q = 10 Myr. Uranus 
begins with a — 14.5 and e = and also undergoes mi- 
gration, to 19.3 AU, on the same timescale. Jupiter and 
Saturn begin with a = 5.2 AU, e = 0.05 and a = 9.6 AU, 
e = 0.05 respectively and do not undergo migration or ec- 
centricity damping. The integra tion includes 24,000 test 
particles, half of which (following lLevison et al.ll2008l ) be- 
gin in the region from 20-29 AU with e = 0.2 and half of 
which begin in the region from 29 - 34 AU with e = 0.15. 
The two example particles shown were among the group 
of particles found in the stable classical region in this in- 
tegration after 4 Gyr and exhibit typical behavior. The 
particles are scattered into the classical region above the 
region of survival (Figure 3; below dotted line in right 
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Fig. 9. — Other giant planets affect the cold classicals indirectly 
by causing apsidal precession of Neptune. Top panel: snapshot 
(1.4 Myr) from an integration including four planets (with initial 
conditions a = 5.2 AU, e = 0.05, i = 1.3° for Jupiter, a = 9.54 
AU, e = 0.06, i = 2.49° for Saturn, and a = 16 AU, e = 0.05, 
i = 0.773° for Uranus). Middle panel: same for integration includ- 
ing just Neptune undergoing apsidal precession with a period of 
1.6 Myr. Bottom panel: same for integration including Neptune 
not undergoing precession. The black line in the top plot is the 
first-order multi-planet secular theory (i.e. without the extra res- 
onant correction terms or higher order terms w e included for 
the just-Neptune theory) (Murray & Dcrmott 2000). The gray line 
on each plot is the analytical expression from the middle panel for 
comparison, computed using single planet secular theory, including 
precession (Equation [4} . The black line in the bottom panel does 
not include precession. Note the scattering in the top two panels 
interior to the 2:1 resonance. 



intersect, r = r^'- 

a(l - e 2 



<xn(1 - e£f) 



1 + e cos(# — w) 1 + eN cos( 



TZ7 N ) 



(10) 



This is important because once the particle is scat- 
tered, it begins to undergo secular oscillations and the 
initial phase of the oscillation, /?o (Equation |4|, depends 
on w — njfj. 

Since in the classical region, most of the particle's orbit 
is outside of Neptune's orbit, the orbits will intersect 
close to the particle's periapse, i.e. the interior part of its 
orbit, so 9 ~ vd. When 9 — vd, Equation (p~0|) simplifies 
to: 



a(l 



a N (l 



1 + cn cos(n7 — TUn) 



(11) 



and thus w maps exactly to the particle's post-scattering 
periapse q. For example, particles scattered to the mini- 
mum q = r p 7v have vd— vd^ = 0, while particles scattered 
to the maximum q ~ r a .N have vd — vd-r = ir. Parti- 



a N (l 



have 



cles scattered to an intermediate 
vd — tUN = i7r/2. 

In Figure 111! we plot the longitude of periapse relative 
that of to Neptune w — versus periapse q of test par- 
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Fig. 10. — Secular evolution can deliver hot objects into the clas- 
sical region. Evolution of two example particles (to p and bottom) 
for the first 6 Myr of an 4 Gyr integration resembling lLevison et all 
( 2008) Run B. Left: path of objects in (a, e) space. The color varies 
from purple (beginning of the integration) to red (6 Myr). The 
dashed line indicates the scattering line q = 35. Right: evolution 
of the particle's eccentricity versus time. The dashed line indicates 
one e-folding time, T eN = 2 Myr, for the damping of Neptune's 
eccentricity. 



tides that were scattered in an integration we performed. 
The w — tun is from the first (3000 yr) timestep after the 
particle's scattering. The integration lasted for 1 Myr 
and included Neptune, with <2n = 30 AU and eN = 0.3, 
and 11,600 test particles evenly-spaced in semimajor axis 
from 29-34 AU, with e = 0.15. The w — otn values are 
well- matched by Equation (fTTj) . 




Fig. 11. — Longitude of periapse relative to Neptune zu — ctn of 
a particle's orbit after being scattered by Neptune into the region 
a > 40 AU maps to the particle's new periapse distance q. The 
particles from t he in tegration are plotted as pluses and the solid 
line is Equation lilt . The dashed lines are the solid line shifted by 
±40°. 



Constraints from classicals 



15 



The secular evolution of the particle after its scattering 
will depend on its semimajor axis, eccentricity, and longi- 
tude of pericenter relative to Neptune's (Equation [2]). In 
particular, from Equation ([2]), it follows that the initial 
rate of change of the particle's eccentricity depends on 
zu — tun- By calculating the particle's total eccentricity 
e = \/h 2 + k 2 , differentiating with respect to time, and 
evaluating the time derivative at t — 0, we find: 

(2) m N n 

e(t = = -eNSin ro — zu^)ab\ L(a) 12) 

i/z m 4 
O 

An analogous expression follows from Equation (|4]). 

Along a scattering line q, the KBOs do not have ran- 
dom zu but zu close to the value dictated by Equa- 
tion (|TT]) . Thus the particles scattered to the minimum 
Q = r p,N, which have zu — zu-r = 0, and maximum 
Q = r a,Nj which have zu — zu-^ — ?r, are turning over 
in their secular evolution cycles (e = 0). However, parti- 
cles scattered to an intermediate q — on(1 — e^j), which 
have zu — zu^ = n/2, will be decreasing in eccentricity at 
the maximum rate in the cycle. 

Moreover, the particle's free eccentricity also depends 
on zu — wn' 

e froo = Coed + (e(0)) 2 - 2e(0)e forced cos(ro - zu N ) (13) 

where e(0) is the KBO's eccentricity at t = 0. 

Thus zu — hjn sets not only the particle's initial phase 
in its secular evolution cycle but also the amplitude of its 
free eccentricity. The particles with the phase to achieve 
the lowest possible total eccentricity (ef rec close to eforccd) 
may not initially be evolving downward in their cycle. 
Therefore, in order to calculate the minimum time for a 
particle to reach the stable region we consider all values 
of zu — TUN- 

Thus particles with zu — zun — and zu — zun — n have 
the smallest and largest free eccentricities respectively, 
while those with zu — zun = ±7r/2 have an intermediate 
value. Unfortunately for particles trying to reach low 
eccentricities, the ones with largest free eccentricity (zu — 
n7N = 7r) are initially going up (i.e. e is increasing). 

From numerical integrations (Figure Hip , it appears 
that the maximum deviation in zu — tun as a function 
of q from Equation (fTTj) is ±40° for conditions relevant 
to Neptune and the classical region KBOs (due to the 
fact that 9 in Equation (fT0|) is not always exactly zu). 
In setting the initial conditions for the secular evolution 
of scattered particles, we employ this mapping between 
zu — ?z7n and the particle's q, including the uncertainty. 

Thus the criteri on fo r delivering the hot classicals de- 
veloped in Section l2~2l requires that: 

CONSTRAINT 2: Neptune's apoapse is large enough 
(^a,N > 34 AU) so that particles are immediately scat- 
tered into the stable region, or Neptune imparts an efo rcc d 
large enough relative to ef rco so that it is possible for par- 
ticles with semimajor axes in the range 42.5-47.5 AU to 
evolve to q > 34 AU in less than r e ^. 

This constraint ensures that it is possible for at least 
some hot objects in the region from 42.5 to 47.5 AU to 
be delivered into the region of longterm stability. 

4.4. Accelerated secular forcing near resonances 

When two bodies are near resonance, the secular eccen- 
tricity forcing happens on a much faster timescale. This 



effect has been recognized as th e ca u se of Saturn and 
Jupiter's fast precession, which iHilll (|1897T ) attributed 
to their period ratios bei ng near 5:2. More recently, 
iMinton fc Malhotral (|2011h recognized that the 2:1 reso- 
nance also contributes to Jupiter and Saturn's fast pre- 
cession and lMalhotra et all ()1989f ) identified the classical 
Uranian satellites' proximity to resonance as the cause of 
their deviation from their predicted ephemerides. 

For the cold objects, the accelerated secular forcing 
near resonance quickly excites the eccentricities of these 
objects (as seen in Figure [7]) , disrupting the confinement 
of the cold population. Near resonance, the correction 
term, e^KBO: to c/kbo is large (Equation [3]) . Thus the 
secular frequency c/kbo is very high. If resonances over- 
lay the cold classical region at early times, Neptune's 
eccentricity would have to damp on unrealistically short 
timescales to fulfill CONSTRAINT 1 (Equation EJ) . As 
shown in Figure [5J the objects near resonance remain 
dynamically disrupted even after Neptune's eccentricity 
damps. Since the cold objects in the region 42.5 AU 
< a < 45 AU are confined to low eccentricities (Section 
12.21) . they cannot have been excited by accelerated sec- 
ular forcing near resonance while Neptune's eccentricity 
was high. 

CONSTRAINT 3: Resonances cannot overlie the re- 
gion 42.5 AU < a < 45 AU while Neptune's eccentric- 
ity is high. This constraint is a special case of CON- 
STRAINT 1 (Equation [9]) and is quantified in Section 
[5X2] 

For the hot objects, accelerated secular forcing near 
resonance can drive down their eccentricities once they 
have been scattered into t he classical region. Figure 
[T2l inspired by Figure 3 of iLevison et al.l (|2008f) . shows 
two example integrations of particles beginning at large 
eccentricities and evolving down to smaller eccentrici- 
ties. The in i tial c onditions (panel i) match those in 
ILevison et al"1(|2008[ ). Figure 3 (their top left panel). The 
snapshots in panels (ii) and (iii) are after 1.4 Myr. The 
integration in panel (ii) includes all four giant planets. 
Neptune has initial conditions <2n = 30 AU and eN = 0.2, 
and the other planets have their modern orbital elements. 
The integration in panel (iii) is the same except with- 
out Uranus. The particles are undergoing secular evo- 
lution, as demonstrated by their vertical paths in (a, e) 
space (a does not change under secular evolution) and 
their eccentricity oscillations (panel iv). If Uranus is 
present near the 2:1 resonance with an eccentric Nep- 
tune, the evolution of the particles is chaotic (panel ii), 
as we will discuss in more detail in Figure 14.51 For ob- 
jects capable of reaching low eccentricities through secu- 
lar evolution (Equation [13]), proximity to resonance sig- 
nificantly decreases the delivery time. Fast secular evolu- 
tion near resonances can also assist in capturing objects 
into resonance, in addition to the previously iden tified 
mechanisms of cha otic capture (ILevison et al.ll2008h and 
smooth migration (lMalhotralll995D . 

4.5. Chaotic sea: no additional constraints 

When Neptune's eccentricity is large, the reso- 
nances are w id ened and potentially overlap in what 
(120081) describe as a "chaotic sea." 
20081 ) argue that this region extends to 
the 2:1 resonance when Neptune's eccentricity is eN > 
0.15. We have investigated the circumstances for chaos 
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Fig. 12. — Objects scattered into the classical region undergo 
re-scattering and secular evolution. Panel (i): Initial conditions of 
the particles. Panel (ii): Positions of the p articles (black circles ) 
after 1.4 Myr (compare to the middle row of Levison ct al. (2008), 
Figure 3) in an integration including all four giant planets. The 
gray is the cumulative region visited by the particles. Panel (iii): 
Same for an integration without Uranus. The colors show the paths 
of a few selected particles. Panel (iv): Eccentricity oscillations for 
particles corresponding to panel (iii). 

and reached several conclusions, which we will state and 
then justify: 

1. Even when widened by Neptune's high eccentric- 
ity, the resonances between the 5:3 and 2:1 do not 
overlap except for particles at high eccentricities. 

2. Variations in Neptune's semimajor axis on 
timescales of order a KBO libration time combine 
with widened resonances to cause a chaotic region. 

3. The existence and extent of the chaotic region de- 
pend on the details of Neptune's interactions with 
Uranus. 

4. The potential for chaos does not impose immedi- 
ate additional constraints beyond those described 
in the previous subsection (Section 14.2.21) . We will 
show that we can still model just Neptune, tak- 
ing into account its potential evolution under the 
influence of the other giant planets. 

4.5.1. Just Neptune: no chaotic sea 

In addition to examining each particle individually, a 
qualitative way to distinguish between: 1) a chaotic sea, 
and 2) particles secularly evolving until they are scat- 
tered, is to plot the individual paths of a collection of 
particles through (a, e) space (Figure [T3l . First we con- 
sider an integration that includes only Neptune (row 1). 



Each starting at e = 0, the particles move straight up- 
ward vertically in (a, e) space until they reach the scat- 
tering line (solid black line, e = 1 — r a ^/a), rather than 
moving horizontally and vertically as they would in a 
chaotic sea. 

In integrations including just Neptune, with or with- 
out apsidal precession, the eccentricity of an initially cold 
particle grows secularly until its orbit crosses Neptune's 
at e > 1 — r a ^/a. For a particle at 42.5 AU, when 
<2n = 30 AU and eN = 0.2, this threshold is e > 0.15. 
After reaching this threshold, the particle then undergoes 
scattering events. Even particles near resonances evolve 
secularly, with the increased secular frequency defined 
by Equation ([5]). These particles appear separated from 
the coherent excitation of the other cold particles (Figure 
[7]) because they undergo secular evolution so quickly and 
because the secular oscillation rate depends quite steeply 
on semimajor axis near resonance. It appears that, for 
particles that begin at low eccentricities, the resonances 
do not overlap even when eN = 0.2. Therefore the con- 
straints from 14.2.21 hold in the just-Neptune case, even 
when the planet is processing. 

4.5.2. Neptune, Jupiter, and Saturn: no chaotic sea 

Adding Jupiter and Saturn in their current configura- 
tions does not generate a chaotic sea (Figure Q21 row 2). 
The particles continue to move upward in (a, e) space 
until they reach the scattering line. The behavior is not 
qualitatively different from the just-Neptune case. 

4.5.3. Neptune and Uranus: chaotic sea 

However, when Uranus is added on its current orbit, a 
chaotic sea appears in the classical region, extending up 
to the 9:5 resonance (Figure fT3l row 3). In the chaotic 
regime, individual particles exhibit chaotic jumps in their 
eccentricity. Some cross from one resonance to another. 
They move horizontally, as well as vertically, in (a, e) 
space. 

Why does adding Uranus create the chaotic sea? Nep- 
tune and Uranus exhibit anti-correlated variations in 
their semimajor axes associated with proximity to their 
2:1 resonance. For the configuration considered here, the 
periodicity of this variation is about 10 4 years, and the 
amplitude for Neptune is about 0.2 AU. This timescale 
is of order the typical libration time of a resonant KBO 
in the classical region. We performed additional in- 
tegrations in which we modified Mercury6 to turn off 
the gravitational interaction between each KBO and any 
planet except Neptune. The behavior was qualitatively 
the same as that shown in Figure [T31 row 1, suggest- 
ing yet again that the other giant planets only affect 
the KBOs indirectly through their impact on Neptune's 
orbit. We performed integrations that include just Nep- 
tune, no other giant planets, with its semimajor axis os- 
cillating with a period of 12,400 yr and amplitude 0.2 
AU. The behavior was the same as in the four planet 
case, and the chaotic sea was present (Figure [T3l row 
4). Evidently this strong periodicity in Neptune's orbital 
variations, which is driven by Uranus and which moves 
the locations of the resonances on timescales of order a 
KBO libration time, causes the chaotic sea. 

We note that in the case in which Uranus is excluded 
(Section 14.5.21 Figure [T3l row 2), Neptune's semimajor 
axis is perturbed on orbital timescales (165 years) by 
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Fig. 13. — Paths of particles in (a, e) space provides a qualitative 
probe of the existence of a chaotic sea; the solid line indicates the 
scattering line. The locations of the resonance centers arc plotted 
as dashed lines; in the case in which Neptune's semimajor axis 
changes (i.e. in rows 3 and 4), the minimum and maximum centers 
are plotted. We integrated 200 massless test particles, each starting 
with e = 0, for 6 Myr under the influence of a subset of planets. 
Neptune, present in each integration, begins with ajvj = 30, ejM = 
0.2. Row 1: the top panel is an integration with just Neptune, 
precessing with a period of 4 Myr. The chaotic sea is not present: 
particles evolve secularly upward in (a, e) space until they reach 
the scattering line. Row 2: the behavior of the particles in this 
integration, which also includes Jupiter and Saturn at their current 
locations, is qualitatively similar: no chaotic sea. Row 3: in this 
integration — which includes Jupiter, Saturn, and Uranus at their 
current locations — the chaotic sea appears, extending up to about 
45 AU, just past the 9:5 resonance. Instead of a straight upward 
path, the particles move from left to right as well. Row 4: the 
chaotic sea is also present in this integration, which includes just 
Neptune but with its semimajor axis oscillating with a period of 
12,400 yr and amplitude 0.2 AU, as it does under the influence of 
Uranus in row 3. Values for a and e are in barycentric coordinates. 



Saturn and Jupiter, with an amplitude of 0.02 AU. This 
small-amplitude perturbation in on, on a timescale two 
orders of magnitude shorter than the resonant timescale 
of KBOs, does not create a chaotic sea. The chaotic sea 
appears to be limited to particular dynamical histories 
in which Neptune and Uranus are strongly interacting 
through the 2:1 resonance. A thorough exploration of 
these histories are beyond the scope of this paper, but 
may provide additional constraints. 

The chaotic sea appears to extend to 45 AU, just past 
the 9:5 resonance. This region between the 3:2 and 9:5 
resonances is already forbidden to overlie today's cold 
classical region because the secular precession rates are 
extremely fast there (Section l4.4l) . so — in the case of the 
cold classicals — the potential for chaos adds no addi- 
tional constraints. We notice that 45 AU coincides with 
the current edge of the cold classical region. In our in- 
terpretation, t his is a coincidence. We also note that 
iLevison et al.l (|2008f) found that the chaotic sea extended 
to the 2:1. Because the chaos depends on the interac- 
tions between Uranus and Neptune, we expect that this 
difference may be due to different initial conditions for 
Uranus. Moreover, their interpretation that the chaotic 
sea extends to the 2:1 is based on the fact that the cu- 
mul ative region visited by the particles extends to the 
2:1 (|Levison et al.l l2TjQ8t, Figure 3). An alternative inter- 
pretation is that chaotic sea extends only to the 9:5 but 
that particles beyond this location are close enough to 
either the 9:5 or the 2:1 resonance to quickly reach low 
eccentricities through secular evolution, which is faster 
in regions near resonances. 

Ou r constraints based on secular evolution (Section 
I4.2.2[) are conservative. The chaotic sea cannot revive 
a region of parameter space which we have excluded, 
but it can rule out additional regions. If oscillations in 
Neptune's semimajor axis caused by interactions with 
Uranus are large enough, the chaotic sea may extend be- 
yond the 9:5 resonance, which would impose additional 
constraints. For example, if the chaotic sea extended 
to the 2:1 resonance, the region of parameter space for 
Neptune with eN > 0.1 between 28-29 AU (which we will 
demonstrate is viable in Section 5.1, Figure fTS"]). in which 
the cold classical region is sandwiched between the 9:5 
and the 2:1 resonance, would no longer preserve the low 
eccentricities of cold classicals. 

Another effect of the oscillations in Neptune's semi- 
major axis is to effectively widen the resonances. This 
effect could potentially cause more KBOs at the edge 
of, but not within, the chaotic sea (for example, KBOs 
just beyond the 2:1 resonance) to experience fast secular 
evolution due to proximity to orbital resonance (Section 
14. 4|) . Thus even more parameter space could be ruled 
out. We do not explicitly take this into account because 
the oscillations in Neptune's semimajor axis depend on 
the particular configuration of Neptune and Uranus. 

Our constraints are a starting point for more extensive 
integrations, which will require careful consideration of 
the interactions between Neptune and Uranus. We note 
that, for illustrative purposes, we have used the current 
semimajor axis and eccentricity of Uranus and current 
semimajor axis of Neptune in these explorations, but the 
existence and extent of the chaotic sea depends on their 
particular orbital configuration — especially their prox- 
imity to the 2:1 resonance — during Neptune's wild days. 
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4.6. Summary 

By modeling the four dynamical processes that result 
from Neptune's high eccentricity — scattering (Section 
14.11) . secular forcing (Section 4.2-4.3), accelerated secular 
forcing near resonances ( Section I4.4[) . and a chaotic sea 
(Section 14. 5[) — we have translated the conservative cri- 
teria imposed by the observed eccentricity distributions 
of the hot and cold classicals (Section 12. 2|) into the fol- 
lowing constraints: 

• Neptune's apoapse must be large enough to deliver 
hot objects to the longterm-stable classical region 
immediately (r a ^ > 34 AU) or ef orC cd must be 
large enough (relative to ef reo ) to evolve the parti- 
cle's e to < 0.3 in less than Neptune's eccentricity 
damping time r eN 

• The final value for the eccentricities 
of planetesimals in the region from 
42.5 AU < a < 45 AU must be less than 
e = 0.1: sin(min( gK Bor eN , tt/2)) >*° bo e N < 0.1 

(or the forced eccentricities must be kept below 
0.1 by fast precession). 

• Resonances cannot overlie the region 42.5 AU < 
a < 45 AU while Neptune's eccentricity is high. 

Thus Neptune's eccentricity must be high in order to 
deliver hot objects to the classical region and yet will 
disrupt the cold objects quickly unless Neptune's high 
eccentricity damps quickly or the planet's orbit apsidally 
precesses quickly. In all cases, mean motion resonances 
with Neptune cannot overlie the region 42.5 AU < a < 45 
AU while Neptune's eccentricity is large. 

5. RESULTS: CONSTRAINTS ON NEPTUNE'S DYNAMICAL 

HISTORY 

Applying the constraints developed in Section [4j we 
place constraints on parameters of Neptune during its 
dynamical history. First, we consider separately which 
parameters of Neptune preserve an in situ cold popula- 
tion (Section I5.1[) and which allow delivery the hot clas- 
sicals (Section 15.21) . In Section 15.31 we combine those 
constraints and determine which parameters of Neptune 
allow the planet to simultaneously preserve the cold clas- 
sicals while delivering the hot classicals. In Section 15.41 
we interpret these parameter constraints in light of Nep- 
tune's full dynamical history. We present example inte- 
grations illustrating the constraints in Section 15.51 

The combined constraints will offer answers to the fol- 
lowing questions about Neptune's dynamical history: 

• Could Neptune have been scattered to a high- 
eccentricity orbit? 

• If so, how quickly did dynamical friction damp 
Neptune's eccentricity? 

• How far did Neptune migrate in the protoplanetary 
disk? 

• If both damping and migration occurred, what 
were their relative timescales? 



5.1. Regions of parameter space that keep cold objects at 
low eccentricities 

In Section 15.1.11 we identify which regions of parame- 
ter space fulfill CONSTRAINT 1 (Equation ^ , preserv- 
ing the cold classicals below e < 0.1 (as summarized in 
Section [476)) . without including the effects of orbital res- 
onances or precession. In Section 15.1.21 we incorporate 
CONSTRAINT 3 th e effects of orbital resonances. Fi- 
nally, in Section [5.1.31 we consider the special case of fast 
precession. 

5.1.1. Constraints on Neptune's eccentricity and damping 

time 

We begin by identifying regions of parameter space 
where, for a given semimajor axis of Neptune on, Nep- 
tune's eccentricity eN is small enough or its damping time 
t £n is short enough to avoid excessively exciting the cold 
classicals. In this subsection, we neglect the effects of 
resonances and assume zero precession of Neptune's or- 
bit, in Figure [T~4] reveals the two two main regions of 
parameter space. 
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Fig. 14. — Constraints on Neptune's parameters to preserve the 
forced eccentricities of KBOs below ef orcc( j = 0.1 in the region from 
42.5 to 45 AU, where aN optune = ajj and eN cptuIle = are spelled 
out for clarity. The white region indicates parameters of Neptune 
that keep ef orce< j < 0.1, no matter how long the damping time. 
The gray regions are contours of the maximum allowed log 10 r eN 
(where r eN is in years), neglecting orbital resonances and assuming 
Neptune's orbit has zero precession, to fulfill CONSTRAINT 1 
(Equation . 



1. In the contoured region (high cn), we plot con- 
tours of maximum eccentricity damping time as a 
function of (aN,eN) that fulfill the criteria set by 
CONSTRAINT 1 (Equation [9]). We calculate the 
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maximum eccentricity damping time by consider- 
ing the equation in CONSTRAINT 1 for a particle 
at 42.5 AU, where the secular evolution is fastest 
(excluding resonances). This map of (on, &n) illus- 
trates constraints on how quickly Neptune's eccen- 
tricity must damp — in order to avoid exciting the 
cold classicals above the level we observe — when 
Neptune occupies a particular region of (on, cn) 
space. 

2. In the white region, Neptune's eccentricity is small 
enough that the forced eccentricity of the cold clas- 
sicals never exceeds ef orcc d = 0.1 and thus the 
damping time can be arbitrarily long. As shown in 
Wolff etaLl [20121 in the regime of slow damping, 
the particle's eccentricity damps to its free eccen- 
tricity which, for a KBO beginning with e = 0, is 
equal to ef orccd . 

We will present several more such plots throughout 
the paper. Note that we plot only eccentricities up to 
eN = 0.4, corresponding to a r p ^ = 18 AU for on = 30 
AU. Above this value, corrections for Neptune's high ec- 
centricity beyond what we already included would be 
necessary, and one also worries about the orbit of Nep- 
tune crossing the orbit of Uranus, which is currently sit- 
uated at 19 AU. However, the constraints we have devel- 
oped could be considered for larger eN- 

The criteria in this subsection hold when the secular 
excitation times are not affected by proximity to reso- 
nance. 

5.1.2. Constraints on Neptune's dynamical history, 
including the effects of resonances 

In regions near orbital resonance with Neptune, KBOs 
undergo significantly faster secular evolution, as demon- 
strated in Section 14.21 Here we incorporate the reso- 
nance correction terms for the secular excitation times, 
in Figure HH analogous to Figure HH we plot contours 
of eccentricity damping time as a function of (on, cn) 
that fulfill the criteria set by CONSTRAINT 1 (Equa- 
tion [9J . We calculate the damping time of Neptune's 
eccentricity, r e for which 80% of initially cold objects in 
the region 42.5 AU < a < 45 AU remain below e = 0.1. 
Note the key difference between the two figures: in Fig- 
ure 1 141 r e varied smoothly through (an, e^) space, but in 
Figure [T51 there are dark regions where the eccentricity 
damping time is substantially reduced due to resonances 
overlying the classical region. 

The damping times are unrealistically short for param- 
eters of Neptune for which resonances overlie the cold 
classical region. Neptune is unlikely to have spent sub- 
stantial time w i th hig h eccentricity in these <zn ranges. 
iFord fc Chiang (|2007f ) nnd that a Neptune -mass planet 
at 20 AU with an eccentricity of 0.3 undergoes eccen- 
tricity damping on r 6N = 0.6 Myr to 1.6 Myr (depend- 
ing on whether Neptune's orbit intersects the planetes- 
imal disk at pericenter/apocenter or at quadrature) if 
the surface density a of planetesimals is 1 g/cc. This is 
roughly the surface density needed to grow plane tesimals 
in the region of Neptune (jKenvon fc Luulll998l ). bmce 
r eN cx o -- 1 , it is unlikely to be less than 0.1 Myr (10 5 yr), 
an order of magnitude faster than the r eN calculated by 
IFord &; Chiand (|2007l) . With resonances incorporated, 
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Fig. 15. — Contours of log 10 T eN (where Te N is in years), from the 
constraint that the small bodies be preserved below e < 0.1 in the 
region from 42.5 to 45 AU, where aNeptune = a N an d eNeptune = 
eN are spelled out for clarity. The white region indicates where 
Neptune's eccentricity is small enough that the forced eccentricity 
of the cold classicals never exceeds e = 0.1. Note that the dark 
regions arc forbidden to Neptune; in these regions, the resonances 
overlie the cold classical region, exciting the cold classicals to high 
eccentricities on timescales less than 10 4 ' 5 years. The constraints 
in this plot correspond to the situation in which at least 80% of 
the particles in the region from 42.5 to 45 AU are confined below 
e < 0.1. 

certain regions cannot satisfy CONSTRAINT 1 (Equa- 
tion [5]) in the zero precession case without un-physically 
low values of r GN . CONSTRAINT 3 is a qualitative state- 
ment of this result. Thus by applying CONSTRAINT 1, 
including resonances, we recover CONSTRAINT 3. 

Incorporating the resonance correction terms, we plot 
the maximum r eN as a function of the KBO's semimajor 
axis for several combinations of aN and eN (Fig [TBI . The 
maximum eccentricity damping time is larger for lower 
eccentricities, as shown for several illustrative values of 
6n in the top panel. The dips where r eN approaches 
correspond to regions of the Kuiper Belt near mean mo- 
tion resonances where the secular excitation time is very 
short. The dips are wider when Neptune's eccentricity 
is higher. When these dips overlie the cold classical re- 
gion from 42.5 AU < a < 45 AU, the damping time 
requirement is un-physically short. In the bottom panel, 
we see how the resonances shift with Neptune's semima- 
jor axis. At aN = 30 AU, several resonances overlie the 
cold classical region. When aN = 28 AU, the resonances 
are shifted interior and the cold classical region is sand- 
wiched between two resonances, the 9:5 and the 2:1. For 
aN = 27.5 AU (not shown), the 2:1 resonance would be 
on top of the cold classicals. When aN = 26 AU, all the 
major resonances are interior to the cold classical region. 
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Fig. 16. — Maximum eccentricity damping time T eN (abbreviated 
Te) to preserve the cold classicals at e < 0.1, as a function of a 
KBO's semimajor axis for illustrative parameters of Neptune. The 
gray region is the cold classical region (42.5 AU < a < 45 AU). 
Top panel: maximum eccentricity damping times for: a>j = 30 
and = 0.4, 0.3, 0.2, and 0.13. Bottom panel: times for a-^ = 
30,28.7,26 for e N = 0.2. 

5.1.3. Special regimes of fast precession 

In this su bsection, we explore the special case of fast 
precession. iBatvgin et al.l ()2011[ ) first suggested that if 
Neptune were to prcccss sufficiently quickly, the cold 
classical population would remain unexcited because fast 
precession lowers the forced eccentricity (Equation[5~J). In 
Figure [TT1 we plot the forced eccentricity as a function 
of the particle's semimajor axis for a range of precession 
rates when on = 30 AU, eN = 0.3. Thus r eN in the 
allowed region can be arbitrarily long if the precession 
period is sufficiently short. A precession period of 0.9 
Myr keeps ef orcc d < 0.1 in the cold classical region for 
these parameters of Neptune (compare to the current gs 
precession period of 2 Myr). However, an even faster 
apsidal precession period for Neptune may be required 
if not only secular evolution but scattering and/or chaos 
excites the cold classicals. 
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Fig. 17. — Forced eccentricity vs semimajor axis for a range of 
precession rates, for eN = 0.3, = 30. The forced eccentricity 
stays below 0.1 in the region 42.5 AU < a < 45 AU for 27r/ro < 0.9 
Myr 



Unlike in the case of fast damping (Section 14.1.11) , in 
the case of fast precession, scattering cold KBOs by Nep- 
tune may impose an important constraint. In the case 



of fast precession (and slow damping) , even though the 
KBO's eventually evolve to ef rec once Neptune's eccen- 
tricity damps, the KBO can reach a maximum value 
of e max = 2ef orce d while eN is still high. Thus the 
KBO is vulnerable to being scattered. For example, 
fora N = 30 AU, e N = 0.3, a KBO at 42.5 AU can be 
scattered if it reaches e > 0.08. Thusthe KBO needs 
eforccd < 5 0.08 = 0.04 to avoid being scattered, requir- 
ing a fast precession period of 0.4 Myr for Neptune. 

Chaos (Figure [T3")) can still excite the cold classicals 
even if their forced eccentricity is low. In this case, res- 
onances cannot overlie the cold classical region if Nep- 
tune's semimajor axis oscillates strongly due to interac- 
tions with Uranus (Section l4.5l) . Therefore, if the chaotic 
sea is present, the dark regions in Figure [TBI are still for- 
bidden. 

Figure [18] demonstrates how fast precession caused by 
Neptune's interaction with the other planets is unable to 
keep the cold objects at low eccentricities in resonance 
regions. We performed two integrations; in each of which 
Neptune has on = 30.06 and eN = 0.2, and 600 test par- 
ticles begin with e — 0. The first integration (top panel) 
includes the other three giant planets, on their current 
orbits; they cause Neptune to undergo apsidal precession 
with a period of 1.2 Myr. In the second integration, Nep- 
tune is forced to precess at this rate without any other 
planets included (see Section [3731 for details on the imple- 
mentation). We plot the maximum eccentricity reached 
by each particle. In the top panel, chaos has dynamically 
disrupted particles that were preserved at low eccentric- 
ities by fast precession in the bottom panel (i.e. in the 
region interior to 45 AU). 
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Fig. 18. — If fast precession is caused by Neptune's interactions 
with Uranus, chaos can disrupt the confinement of the cold clas- 
sicals. Top panel: maximum eccentricity reached by each particle 
during a 1.6 Myr integration including all four giant planets. De- 
spite fast precession, objects interior to the 2:1 resonance are ex- 
cited to high eccentricities. Bottom panel: maximum eccentricity 
when just Neptune precesses at the same rate (period of 1.2 Myr) 
as in the four-planet case. 

Unfortunately, the most obvious configuration that 
causes fast precession — Neptune and Uranus near their 
2:1 resonance — also causes large oscillations in Nep- 
tune's semimajor axis. The larger the oscillations in 
Neptune's semimajor ax i s, th e wider the chaotic sea. 
Although Bat ygin et al.l (|2011l) suggest that a massive 
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planetesimal disk could contribute to fast precession, this 
precession rate would be roughly mdi3k ra^, where tin is 

Neptune's orbital frequency. Thus a disk mass of roughly 
60 Earth masses would produce a precession period of 0.9 
Myr. Therefore, it is unlikely that the planetesimal disk 
alone could produce sufficiently fast precession. When 
Neptune is in the region from 28-29 AU - and the clas- 
sical region is sandwiched between the 9:5 and 2:1 reso- 
nance - oscillations in Neptune's semimajor axis due to 
resonant interactions with Uranus could cause a chaotic 
sea that extends into this "sandwich" region. Thus fast 
precession may fail for parameters for which fast damp- 
ing has the possibility of working. 

We will present a contour map for Neptune's precession 
rate, analogous to Figure in Section |5"U1 

5.2. Constraints on transporting the hot objects to the 
classical region 

Now we consider which parameters of Neptune will al- 
low the transport of the hot objects from the inner disk 
into the classical region. We consider the full range of 
(a, e) into which Neptune can scatter objects from the 
inner disk and the resulting secular oscillations of the 
objects once they reach the classical region. We calcu- 
late a minimum eccentricity damping time as a function 
of (on, cn) by requiring that particles can reach the sta- 
bility region defined by q > 34 AU (CONSTRAINT 2). 
The "minimum time" (contoured) is the time at which it 
is possible for objects in 50% of semimajor axes intervals, 
Aa = 0.063 AU, to reach this stable region. 

Figure [THl a map of constraints on Neptune's param- 
eters that deliver the hot objects, illustrates that there 
are three outcomes for the transport of hot objects: 

1. White region: If r a ,N is sufficiently large, some ob- 
jects are scattered into the stable region. This cri- 
terion can be quantified as r a n = <Xn(1 + 6n) > 34 
AU. 

2. Black region with red stripes: If r a .N is too small, 
the KBOs are scattered into a region that is not sta- 
ble over 4 Gyr and their forced eccentricity is too 
small (relative to ef rcc ) to allow them to secularly 
evolve down into the stable region. This second 
outcome can be understood from Equation |(3J). If 
the forced eccentricity is small, the free eccentricity 
will be close to the particle's initial (large) eccen- 
tricity after scattering. Thus the free and forced 
eccentricity vectors cannot destructively cancel to 
achieve a low total eccentricity. 

3. Middle region with gray-scale timescale contours: 
For intermediate values of r Q ,N , the particles are 
scattered above the stable region (i.e. to eccen- 
tricities too large for stability) but can secularly 
evolve down into the stable region if Neptune's ec- 
centricity damping timescale is long enough. The 
time contours refer to the minimum time (defined 
above) for particles to secularly evolve to the sta- 
ble region. We note that, because the analytical 
model is only first order in the KBO's eccentric- 
ity, the timescales in the intermediate region may 
be in error by up to a factor of two (see Appendix 
14.2.21 for discussion). However, the particles most 
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Fig. 19. — Parameter space of Neptune consistent with delivering 
the hot objects from the inner disk to the classical region, where 
a Neptune = a N an d ^Neptune = e N are spelled out for clarity. The 
white region in the upper right is where particles are immediately 
scattered into the stable region, q < 34 AU. The bottom left region 
(black with red stripes) is where particles can never reach the sta- 
ble region because their forced eccentricity is too small relative to 
their free eccentricity. In the contoured region (middle), particles 
can secularly evolve into the stable region after being scattered if 
Neptune's eccentricity remains high for long enough. The contours 
represent the minimum time (see the text) for any scattered KBOs 
to evolve into the stable region. The pluses mark the conditions of 
the top panel and the bottom panel of Figure [6] which illustrates 
the region into which Neptune can scatter objects for two sets of 
illustrative parameters. 

relevant for this scenario — those rapidly declining 
in eccentricity — are matched with substantially 
smaller error. We also note that this regime (Out- 
come 3) is a narrow region of parameter space and 
that most parameters for Neptune fall within one 
of the two regimes above. 

The minimum values for ((ZNjCn) in the intermediate 
range of r a ,N has a messy analytical solution. The mini- 
mum eccentricity a particle can achieve is: 

e m in = ( e frco + e forccd — 2ef rcc eforccd) (14) 

Combining with Equation (jllj) and (|13j) . for a given aN 
one can solve for the cutoff eN, below which particles 
cannot evolve into the long-term stable region. 

5.2.1. Effects of resonances and chaos 

As shown in Section l5.1.2l the secular forcing near res- 
onance is much faster than in regions outside of reso- 
nance. Here we address the effects of resonances for the 
three outcomes we discussed above: 

1. If r Qi N is sufficiently large, some objects are scat- 
tered directly into the stable region. CON- 
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STRAINT 2 (q > 34 AU) can always be met in 
this region, regardless of secular evolution, so we 
consider this region to always work for the hot 
classicals. Mean motion resonances can make this 
already-allowed region even better, allowing even 
more objects to reach the stable region. 

2. If r 0) N is small, objects are scattered to eccentrici- 
ties that are too large to allow the objects to reach 
the stable region, even through fast secular evolu- 
tion near resonance. Although fast secular evolu- 
tion near resonance lets objects reach their mini- 
mum eccentricity (Equation I14[) . unfortunately all 
the objects have minimum eccentricities too large 
for long-term stability. 

3. In the intermediate regime (the small contoured 
strip of parameter space in Figure Ql|] that does not 
deliver the hot classicals "directly" or "never"), ob- 
jects are not scattered directly to the stable region 
but can reach the stable region through secular evo- 
lution, including fast secular evolution near reso- 
nances. Without resonances, the secular evolution 
timescale at 47.5 AU (the largest KBO semimajor 
axis) sets the minimum timescale over which Nep- 
tune needs to maintain its high eccentricity. How- 
ever, because of fast secular evolution near reso- 
nance, the timescale is actually set by the largest 
KBO semimajor axis unaffected by proximity to 
resonance. Thus the minimum timescale is shorter 
(i.e. Neptune can damp more quickly) than if we 
neglected the effects of resonances. In practice, the 
secular evolution time is n ot a steep functio n of the 
particle's semimajor axis ([Wolff et al.112012 . Figure 
4). We created another version (not shown) of Fig- 
urc[T5]without including the effects of mean motion 
resonances. We required objects at 47.5 AU to be 
able to reach e < 0.3, because this would ensure 
that objects interior of 47.5 would also reach the 
stable region. The results were qualitatively the 
same, with some slightly longer minimum times. 

Now we consider whether objects that are actually in 
resonance, as opposed to being near resonance, could 
be more easily delivered to the classical region than our 
constraints would imply. Although we did not explicitly 
take into account objects actually librating in resonance, 
we argue that, as a by-product, our constraints include 
this effect. When Neptune's eccentricity is low or its 
semimajor axis places no strong resonances in the clas- 
sical region, objects actually-in-resonance with Neptune 
occupy a very small portion of phase space. In the oppo- 
site situation, in which Neptune's eccentricity is high and 
strong resonances are present in the classical region, the 
parameters for Neptune are already in Outcome 1 above 
("directly" scattered). For exam ple, the initial condi- 
tions used to create the Figure 3 of iLevison etH] (2008) 
are part of this region. The constraints we will place in 
this region will come only from the cold objects. 

Potentially the 2:1 resonance could "sweep" the clas- 
sical region, rapidly evolving transported objects into 
the stable region, as Neptune migrates. Moreover, res- 
onances sweeping the region could deliver hot objects 
to the stable region via orbital chaos. However, neither 
of these scenarios would be consistent with preserving 



the cold classicals. Because the hot objects span the 
whole region from 42.5 AU < a < 47.5 (Figure 3, right 
panel), the 2:1 resonance would have to sweep all the 
way through this range, exciting the cold classicals (as 
would the 9:5 resonance as it passes 42.5 AU). Moreover, 
if interactions of Neptune and Uranus created a chaotic 
sea (Section |4.5|) . the cold classicals would be dynami- 
cally disrupted. Therefore, the location of resonances do 
not provide additional constraints or possibilities for the 
transport of hot classicals. 

5.2.2. Effects of precession 

Precession has several effects on delivery of the hot 
classicals: 

1. Before the eccentricity of Neptune damps, preces- 
sion allows Neptune to scatter more objects from 
the inner disk. If a KBO in the region near Neptune 
is on an eccentric orbit, the orbit of Neptune and 
the small body will not necessary cross, depending 
on the relative positions of their periapses. But 
if Neptune's periapse quickly precesses, Neptune's 
orbit may come to intersect additional KBO orbits. 

2. As an object secularly evolves after being trans- 
ported to the classical region, its longitude of pe- 
riapse varies. Hence, precession of Neptune's orbit 
changes the likelihood that the object will be scat- 
tered again before it reaches the long-term survival 
region. We note that the likelihood of scattering 
depends on the relative precession rates of Neptune 
and the object. 

3. If the precession rate is fast enough, it decreases 
the forced eccentricity. The magnitude of the par- 
ticle's free eccentricity is given by Equation (TIB")) . 
For particles with quickly decreasing eccentricity 
(m — tz7n = tt/2) (Equation I12p, a reduction in 
eforced decreases the amplitude of the free eccen- 
tricity. Thus fewer particles scattered above the 
stable region are ever able to reach it. Certain 
(on ; eN ) that allowed particles to evolve into the 
stable region, when Neptune was not quickly pro- 
cessing, will no longer work. 

4. The precession of the forced eccentricity allows par- 
ticles that are evolving down into the stable re- 
gion via secular oscillation to reach the region more 
quickly. This effect changes the timescales for ob- 
jects evolving to the stable region. 

Thus precession can either help or hurt the transport 
of hot objects into the classical region, depending on the 
particular combination of parameters. 

5.2.3. Scattering efficiency of the hot objects 

The number of hot objects that reach the stable region 
depends on the surface density of planetesimals in the in- 
ner disk, the eccentricity damping timescale of Neptune, 
and Neptune's precession rate. A faster precession rate 
evidently delivers more objects to the classical region, 
but objects are also more likely to exit the classical re- 
gion via subsequent frequent scatterings, resulting in a 
steady state flux. Perturbations from other planets could 
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increase the efficiency by increasing Neptune's precession 
rate and causing additional perturbations to the inner 
disk. We leave a more detailed exploration of the scat- 
tering efficiency for future work. 

However we contrast the inefficient process of trans- 
porting the hot objects into the classical region to the 
highly efficient preservation of the cold objects, i.e. all 
the cold objects stay in the classical region but few hot 
objects are transported there. Such a discrepancy might 
be possible if the inner disk is much denser than the outer 
disk. A dense inner disk is consistent with the short 
eccentricity damping timescales we are finding (Section 
15. ip . because the dampi ng rate scales with t he surface 
density of planetesimals (|Ford &: Chi ang 20071) . A dense 
inner disk and rarefied outer disk may explain why Nep- 
tune ceased its pla netesimal-driven migration when it 
reached 30 AU (see iLevison & Morbidellil l200l for the 
suggestion that the planetesimal disk was truncated at 
30 AU). However, the low number of cold objects poses a 
problem for their formation, as we will discuss in Section 

m 

5.3. Combined constraints from both hot and cold 
objects 

We placed constraints on parameters of Neptune that 
preserve the con finem ent of cold objects to low eccen- 
tricities (Section 15. ip or allow the transportation of the 
hot objects from the inner disk into the classical region 
(Section I5.2[) . . These constraints are useful separately, 
and they confirm in situ formation as a feasible origin 
for the cold objects and transport from the inner disk as 
a feasible origin for the hot objects. In this Section, we 
investigate which parameters permit the combination of 
these two origins, producing both a cold and a hot pop- 
ulation. This may be possible if there is overlap between 
the parameter space that preserves the cold classicals 
(Section 15. ip and the parameter space that transports 
the hot classicals (Section IQ)) . 

In Figure 1201 we combine the constraints from the hot 
and cold classicals. First, as a function (<XN,eN) we plot 
the contours for the maximum eccentricity damping time 
necessary to keep the cold classicals at e < 0.1 in the re- 
gion from 42.5 to 45 AU. In the white region, Neptune's 
eccentricity is low enough such that the final eccentrici- 
ties of the cold objects will be below 0.1, no matter how 
long the damping timescale. Then we overplot red, diag- 
onal stripes in the region where Neptune cannot deliver 
the hot classicals, neither by direct scattering nor post- 
scattering secular evolution. In creating this region, we 
took into account the sliver of parameter space (Figure 
[6l contoured region) for which Neptune cannot scatter 
objects directly into the region q > 34 AU but for which 
objects can evolve into the stable region through secu- 
lar evolution. We compare the time required for the hot 
objects to reach the long-term stable region to the maxi- 
mum eccentricity damping time to preserve the cold clas- 
sicals. In most cases, Neptune's eccentricity must damp 
before any hot objects reach the long-term stable region. 
From all these considerations, there is only a small re- 
gion of parameter space where Neptune can deliver the 
hot classicals (no red, diagonal stripes) while fast secu- 
lar evolution near resonances does not quickly excite the 
cold classicals (light regions). 

The allowed region in (on, en) space, in which the hot 
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Fig. 20. — Combined constraints from the hot and cold classi- 
cals when the cold classicals are preserved by Neptune damping 
quickly, where aj| eptuIle = and ejs[ e p tune = are spelled out 
for clarity. The contours refer to the maximum eccentricity damp- 
ing time log 10 T eN (where T epJ is in years) for Neptune in this region 
to avoid excessively exciting the cold classicals. In the white re- 
gion, the forced eccentricity imparted by Neptune in the region 
from 42.5 to 45 AU is below 0.1. In the red, diagonally striped 
region, Neptune cannot transport the hot objects to the classical 
region, defined as q > 34 AU from 42 to 47.5 AU. 

objects can be transported without disrupting the cold 
population, is bounded by a strip extending from (23 AU, 
0.4) to (30 AU, 0.10) on the left and (27.5 AU, 0.4) to 
(30 AU, 0.2) on the right. To the top right of this strip, 
the series of resonances extending from the 5:2 to the 
9:5, widened by Neptune's high eccentricity, are on top 
of the cold classical region. Within this bounding strip 
is a forbidden region near on = 27.5 AU, at which the 
2:1 resonance overlies the cold classical region. 

Based on the considerations above, lower limits can be 
placed on Neptune's migration timescale while Neptune's 
eccentricity is high. During Neptune's period of high ec- 
centricity, Neptune's migration timescale must be slow 
enough to keep Neptune in a region that will not quickly 
excite the cold classicals. Neptune should not spend sub- 
stantial time with high eccentricity near on = 27.5 AU or 
on = 30 AU, lest resonances disrupt the classical region. 
Thus Neptune is constrained to migrate no more than 
a few AU during Neptune's eccentricity damping time. 
Because of the discrete ranges of consistent semimajor 
axes, when resonances are accounted for, a "damp first 
and then migrate" scenario is consistent with preserving 
the cold classicals, while a "migrate first and then damp" 
scenario is not. 

In Figure [2U we plot analogous constraints for the 
special case of fast precession. The contours represent 
the precession rate of Neptune necessary to keep ef orC ed 
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sufficiently low (Equation [7J . In all cases, the KBOs' 
forced eccentricities must be below ef orce d < 0.1 (CON- 
STRAINT 1). We impose an additional constraint to 
avoid the scattering of cold particles as they reach their 
maximum eccentricities of e max = 2ef orcG( j- We require 
e max must be below the scattering line at the inner 
edge of the cold classical region: therefore ef orcc d < 
|(1 — g Qj N/42.5). For large values of Neptune's apoapse 
<7 a ,N, this constraint is stricter than ef orcc d < 0.1. This 
constraint may, in fact, be too strict because if a cold 
classical KBO is scattered by Neptune, it is more likely 
to scattered out of the classical region altogether than to 
end up between 42.5 AU < a < 45 AU at an eccentric- 
ity too large to be consistent with the observations but 
small enough to survive over 4 Gyr. We leave a detailed 
investigation into the role of scattering in the case of fast 
precession for future work. The red, diagonally striped 
region is where hot objects cannot reach the long-term 
stable region. 

Next we consider the constraints for delivering the hot 
classicals. For each (on^n); we calculate whether ob- 
jects can be directly scattered into the region q > 34 AU. 
If not, we assume Neptune imparts a forced eccentricity 
small enough to be consistent with the constraints above 
for keeping the cold classicals unexcited (i.e., eforccd = 
0.1, or, if smaller, efo rco d < |(1 — <7<i,n/42.5)) and deter- 
mine whether the hot object can reach the stable region 
through secular evolution, using Equation ([TJJ . If Nep- 
tune cannot directly scatter objects into the stable region 
and if objects secularly evolving cannot reach q > 34 AU 
in more than 50% of semimajor axis intervals, (<in, Cn) 
is not consistent with delivering the hot objects. In 
the cases where fast precession is necessary to keep the 
cold classicals confined to low eccentricities, hot classi- 
cals scattered into the region have the same low forced 
eccentricities; therefore, they do not experience signifi- 
cant secular evolution down into the classical region but 
remain at their post-scattering eccentricities. 

Finally, in blue vertical stripes, we overplot ranges of 
on for which the centers of one or more resonances of 
fourth order or lower lie in the cold classical region from 
42.5 to 45 AU. These <zn are not necessary forbidden, 
but chaos can occur if Neptune's eccentricity is high and 
interactions with Uranus cause oscillations in a^. 

Though we are applying constraints that are lenient 
and conservative in ruling out regions of parameter space, 
a substantial fraction of parameter space is ruled out. 
The following scenarios have not been ruled out and may 
work, depending on the details of the interactions be- 
tween Uranus and Neptune. In both of these scenarios, 
resonances cannot overlie the region of 42.5 AU < a < 45 
AU. 

• Short T e and large r p j\r: An eccentric Neptune 
transports objects into the classical region until its 
eccentricity damps, which occurs quickly enough 
that the cold classicals are not excited (Figure [5D|) . 

• Fast ton and large r p jy: An eccentric, quickly pro- 
cessing Neptune transports objects into the classi- 
cal region and its quick apsidal precession keeps the 
forced eccentricity of the particles low, preserving 
the cold classicals (Figure l2"Tj). 

In principle, an intermediate r p ^ could deliver hot 
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Fig. 21. — Combined constraints from the hot and cold clas- 
sicals when the cold classicals are preserved by Neptune pre- 
cessing quickly, where a Ncptunc = a N and e Ncptunc = e N are 
spelled out for clarity. The contours refer to the maximum pre- 
cession period r ro = 27r/ri7iM for Neptune to avoid excessively ex- 
citing the cold classicals. Neptune must precess fast enough to 
keep ef orce< j < 0.1 and, when Neptune intrudes into the region, 
e forccd < §(1 ~ 9h,n/42.5). In the white region, the forced eccen- 
tricity imparted by Neptune in the region from 42.5 to 45 AU is 
below 0.1 (and below i(l — q a ^/A2.b), even when Neptune does 
not precess). In the red, diagonally-striped region, Neptune cannot 
transport the hot objects into the stable classical region, defined 
as q > 34 AU from 42.5 to 47.5 AU. The blue, vertically striped 
regions denote a^j for which the center of a resonance, of fourth 
order or lower, lies in the region from 42.5 to 45 AU. These regions 
are not necessary forbidden (see the text). 

objects to high eccentricities (e > 0.3) in the classical 
region but allow them to evolve to lower eccentricities 
(e < 0.3) before the cold objects are excited. In prac- 
tice, we found the timescales are not compatible: if the 
eccentricity damps quickly enough to preserve the cold 
classicals from secular excitation, there is not time for 
hot objects to secularly evolve down into the Belt. How- 
ever, for parameters of Neptune for which hot objects are 
scattered directly in the stable region, secular evolution 
can allow the objects to reach even lower eccentricities, 
especially the objects undergoing fast secular evolution 
near resonances. 

Therefore, in practice, the consistent dynamical his- 
tories are ones in which Neptune has a large enough 
apoapse to transport hot objects immediately into the 
stable region, with its eccentricity damping quickly 
enough or precessing quickly enough so that the cold 
classicals remain at low eccentricities consistent with the 
observations. 

5.4. Interpretation of constraints in light of Neptune's 
full dynamical history 
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The goal of this paper is to determine which param- 
eters for Neptune allow the planet to deliver the hot 
classicals from the inner disk into the classical region 
without dynamically disrupting the in situ cold classi- 
cals. Throughout its dynamical history, Neptune must 
satisfy the constraints presented here to avoid excessively 
exciting the cold classicals (Section 15. ip . At some point, 
Neptune must also spend time in a region of parameter 
space where it can also deliver the hot classicals (Section 
15. 3|) . We clarify that Neptune did not necessary form 
at the location where hot classical delivery takes place 
or arrive there after undergoing a single, instantaneous 
scattering. Before and after hot classical delivery, Nep- 
tune can potentially spend time in any region of param- 
eter space as long it obeys our constraints against not 
excessively exciting the cold classicals. The constraints 
developed here, which can serve as a "road map" for 
Neptune's path through parameter space, hold for real- 
istic models that include multiple scatterings and have a 
straightforward interpretation. 

In the context of the Nice model, Neptune may have 
undergone a series of scatterings, spending time at a va- 
riety of spots in (cin, 6n) space. In each of these spots, 
Neptune must obey the constraints we place to avoid 
disrupting the cold classicals. Perhaps the scattering oc- 
curs quickly compared to the excitation time for the cold 
classicals. If not, Neptune could pre-excite the cold clas- 
sicals (but only to below the observational limit) before 
it reaches a region where it can deliver the hot classicals. 
In imposing our conservative constraints, we are assum- 
ing that the cold classicals begin with e = 0, but if they 
are pre-excited, the constraints will be stricter. 

Another possibility is that Neptune may spend a long 
time at a location where it creates a very small region 
of stability in the classical region, clearing out most of 
the cold classicals. This is a potential solution to the 
mass efficiency problem, which we will discuss in Section 
However, as shown in Section [3.1.11 Neptune cannot 
deliver the hot classicals in this regime. Therefore, after 
this period has ended, when Neptune is delivering the hot 
classicals, the constraints we will place on not exciting 
the cold classicals will hold. 

The scattering(s) Neptune undergoes are quick changes 
in its orbit. After its period of high eccentricity ends - 
or during a temporary period of low eccentricity - Nep- 
tune can undergo slow evolution, including slow migra- 
tion and slow damping of its (now small) eccentricity. 
The KBOs will maintain their free eccentricities through- 
out this slow evolution. If Neptune's eccentricity were 
excited very gradually, on a timescale much longer than 
the secular evolution times of the cold classicals, the cold 
classicals would keep their initial low free eccentricities. 
However, we expect the excitation of Neptune's eccen- 
tricity via scatterings to ta ke place on a timescale shorter 
than millions of years (e.g. pThommc s et al.lll999D . 

We have ruled out much of parameter space with the 
constraints that Neptune cannot excite the cold objects 
above e = 0.1 and must be able to deliver at least a 
few hot objects to q > 34 AU in the region from 42-47.5 
AU. We note that we are "ruling out" parts of parame- 
ter space where Neptune cannot deliver the hot objects 
without disrupting the cold, not "ruling out" that the 
planet can ever spend time there (see above). In Fig- 
ure [22l we look for greater consistency with the obser- 



vations (Figure UJ): a forced eccentricity less than 0.075 
for the cold objects and Oat > 36 AU (meaning that a 
hot object could be scattered to an eccentricity as low 
as 0.24 at 47.5 AU). The parameter space shrinks where 
Neptune can deliver the hot objects without disrupting 
the cold (i.e. light gray regions with no red, diagonal 
constraints). Over-plotted on this figure is an example 
(arrows) of Neptune's path through parameter space. In 
this example, Neptune first undergoes multiple scatter- 
ings, spending a short enough time at each (onjCn) to 
avoid exciting the cold classicals. Then it reaches 28 AU 
with an eccentricity of 0.3; here it delivers the hot clas- 
sicals, as its eccentricity damps quickly enough to avoid 
exciting the cold classicals. Then it migrates at a low- 
eccentricity to its current location. 

On this figure, we overplot some parameters as sym- 
bols. For the triangles and circles, we will show example 
integrations in Section 15.51 The pluses are parameters 
taken from the literature. The plus at on = 20, eN = 0.0 2 
marks the initial condition of lHahn "fc Malhotral (|2005|) . 
from which Neptune undergoes migration to its current 
location. Neptune remains in a region that is white 
(never excites the cold classicals) but also red, diago- 
nal striped (does not deliver the hot classicals) ; thus this 
simulation maintained the low eccentricities of the cold 
classicals but did not deliver the hot classicals. The plus 
at on = 28.9, e N = .3 indicates the initial condition for 
iLevison et alJ (|2008[ ). Run B. In this part of parameter 
space, Neptune can deliver the hot classicals, which were 
indeed produced by this simulation. However, Neptune's 
eccentricity damped on a timescale of r e = 2 Myr, too 
slow to avoid exciting an in situ cold population. 

5.5. Example integrations illustrating constraints 

Example integrations illustrating the constraints we 
have derived are shown in Figure [22] and [M] Each in- 
tegration lasts for 1.6 Myr unless otherwise noted (see 
Sect ion 13731 for a general description of the integrations). 
In addition to the 600 test particles in the region from 
40-60 AU, we add 59600 test particles in the region from 
18-38 AU, each with initial e = i = 0, representing the 
inner planetesimal disk from which the hot classicals are 
scattered. 

The parameters for the four integrations shown in Fig- 
ure [23] are plotted in the parameter space map in Fig- 
ure [22] as circles and correspond to regions of parameter 
space where Neptune cannot both deliver the hot objects 
and keep the cold objects at low eccentricities. In row 1, 
Neptune is at on = 24, cn = 0.02 and undergoes no ec- 
centricity damping. The cold classicals remain confined 
to low eccentricities, but the hot classicals are delivered 
to high eccentricities and cannot secularly evolve down to 
lower eccentricities. For the same reason (i.e. that Nep- 
tune's apoapse is too small), the hot classicals are also 
not delivered in row 2 (<in = 22, eN = 0.35). The damp- 
ing timescale in this integration is longer (2 Myr) than 
half a secular evolution time, so the eccentricities of the 
cold objects converge to the forced eccentricities, which 
are above 0.1. This integration lasts 6 Myr. In the third 
row (<zn = 27.5, eN = 0.35, r e = 0.2 Myr) and fourth row 
(on = 30.06, eN = 0.35, r e = 0.2 Myr), the hot classicals 
are delivered but the cold classicals are excited by fast 
secular evolution near resonances. 

The parameters for the two integrations shown in Fig- 



2G 



Dawson and Murray-Clay 




18 20 



no 
hots 



colds 
okay 



I ' 1 '■ I 
22 24 26 
aNeptune (AU) 



28 30 



<4.5 5 5.5 6 6.5 > 7 
colds okay for log 10 x eN 



Fig. 22. — Combined constraints from the hot and cold classi- 
cal, where aNeptune = a N an d e Neptune = e N are spelled out for 
clarity. The contours refer to the maximum eccentricity damping 
time log 10 T eN (where r eN is in years) for Neptune in this region to 
avoid excessively exciti ng the cold classicals. The constraints are 
stricter than in Figure 1201 In the white region, the forced eccen- 
tricity imparted by Neptune in the region from 42.5 to 45 AU is 
below 0.075. In the red, diagonal striped region, Neptune cannot 
transport the hot objects to the classical region, defined as q > 36 
from 42.5 to 47.5 AU. The arrows are a schematic illustration of an 
example of Neptune's path through parameter space as it under- 
goes multiple scatterings, having its eccentricity damped, or being 
re-scattered, on the contoured timescale to avoid excessively excit- 
ing the cold classicals. T he pluses mar k the initial condition s from 
IHahn fc Malhotral H2005I 1 (bottom) and lLevison et al.l H2008f ) (top). 
The ci rcle s are the parameters for example integrations shown in 
Figure [23] and the triangles for example integrations shown in Fig- 
ure [24] 

ure [24] are plotted in the parameter space map in Figure 
l22l as triangles and correspond to a set of parameters 
in each of the two viable regions of parameter space. In 
both cases, Neptune can deliver the hot classicals. More- 
over, since it is in a region of parameter space where no 
resonances overlie the region of 42.5 AU < a < 45 AU, 
and since its eccentricity damps quickly enough to obey 
our constraints, the cold classicals remain at e < 0.1. In 
the top panel, we plot the observed objects from Figure 
3 for comparison. In the middle panel, the 2:1 resonance 
is interior to the cold classical region. In the bottom 
panel, the cold classicals are sandwiched between two re- 
gions where the cold classicals are excited by fast secular 
evolution near resonance. 

6. DISCUSSION 

Through an exploration of parameter space (Section[5]) 
combined with conservative criteria from the observed ec- 
centricity distributions of classical KBOs (Section I2.2[) . 
we reach the conclusion that most parameters for Nep- 
tune are inconsistent with both delivering a hot popu- 
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Fig. 23. — Examples of four nu mer ical integrations that violate 
the constraints established (Figure [20ll . The blue and red triangles 
are the positions of particles at the end of the integration. The 
red triangles are cold objects which began in the region from 40 
to 50 AU with e = 0. The blue triangles are hot objects which 
began in the inner disk interior to 38 AU. The small triangles have 
eccentricities above the region of long-term stability (dashed line 
in Figure [3J and thus are not expected to survive over 4 Gyr. The 
yellow line indicates e = 0.1 and the shaded gray region is 42.5 AU 
< a <45 AU, where the cold classicals are observed (Figure [3} to 
be confined to low eccentricities. The parameters for Neptune in 
each integration are: on = 24, ejM = 0.02, no eccentricity damping 
(top row), ajj = 22, = 0.35, T e =2 Myr (row 2), apj = 27.5, ejj = 
0.35, T e = 0.2 Myr (row 3), a N = 30, e N = 0.35, r e = 0.2 Myr (row 
4). The snapshots are at times 1.6 Myr, 6 Myr, 1.6 Myr, and 1.6 
Myr. 



lation from the inner disk and preserving a cold popu- 
lation formed in the outer disk. We have explored the 
full parameter space for a generalized model — in which 
Neptune undergoes some combination of high eccentric- 
ity, migration, and/or precession and delivers the hot ob- 
jects on top of the cold — that encompasses the previous 
"chaotic capture" and "extensive migration" models and 
accounts for the different physical properties of the hot 
and cold classicals. We have found that the generalized 
model remains viable only in two restricted regions of 
parameter space: Neptune is scattered onto an eccen- 
tric orbit with a semimajor axis in one of two ranges, 
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Fig. 24. — Examp les of two integrations that obey the constraints 
established (Figure |20| |. The top row shows the observed objects 
from Figure [3] for comparison. In rows 2 and 3, blue and red 
triangle are the positions of particles at the end of the integration. 
The red triangles are cold objects which began in the region from 
40 to 50 AU with e = 0. The blue triangles are hot objects which 
began in the inner disk interior to 38 AU. The small triangles have 
eccentricities above the region of long-term stability (dashed line 
in Figure(3j and thus are not expected to survive over 4 Gyr. The 
yellow line indicates e = 0.1 and the shaded gray region is 42.5 AU 
< a <45 AU, where the cold classicals are observed (Figure [3j to 
be confined to low eccentricities. The parameters for Neptune in 
each integration are: ajj = 26, e^j = 0.38, r e = 0.33 Myr (second 
row) and (In = 28.5, ejM = 0.33, r e = 0.25 Myr (third row). 

24-27 AU (so that the 2:1 resonance is interior to 42.5 
AU) or 28-29 AU (so that the region from 42.5 to 45 AU 
is sandwiched between the 9:5 and the 2:1 resonance). 
Although <2n = 30 AU appears feasible on the param- 
eter space plots for eN ~ 0.15, in this region "strips" 
are excited by overlying resonances and the hot objects 
are far from low enough. Neptune scatters the hot ob- 
jects from the inner disk into the stable classical region 
where we observe them. Because Neptune's eccentricity 
damps or the planet's orbit precesses quickly, Neptune 
does this without exciting the cold objects above their 
observed eccentricities. Because Neptune is confined to 
one of these two particular regions, mean motion reso- 
nances — which would quickly excite the cold classicals 
through accelerated secular forcing and/or chaos — do 
not, while Neptune's eccentricity is high, overlie the re- 
gion where we observe the cold classicals confined at low 
eccentricities today. Most likely, once Neptune's eccen- 
tricity damps, it migrates on a circular orbit to its current 
location, a migration distance of Aon = 1-6 AU. 

Our constraints should be interpreted in light of Nep- 
tune's full dynamical history, which may include multi- 
ple scatterings of Neptune and/or excitation/sculpting 
of the cold classical region before Neptune delivers the 



hot classicals. Throughout its path through parameter 
space, Neptune must obey our constraints on not exces- 
sively exciting the cold classicals. Whatever the prior 
early evolution, Neptune must eventually spend time in 
a region of parameter space where it can deliver the hot 
classicals while its eccentricity damps or precesses quickly 
enough — or Neptune is re-scattered quickly enough- to 
avoid exciting the cold classicals. Then it can proceed to 
its current location via additional scatterings or migra- 
tion, maintaining an eccentricity low enough to continue 
to avoid exciting the cold objects. 

The viable regions of parameter space are qualita- 
tively and quantitatively different from the previous 
models that did not produce the o bserved ec c entric - 
ity distributions. Comp ared to the iMalhotral |l995); 
lHahn fc Malhotral (|2005| ) model, Neptune undergoes a 
period of high eccentricity and migrates a shorter dis- 
tance (< 6 AU, as opposed to 7-10 AU). The m ost sig- 
nificant differences from the iLevison et al.l (|2008f) model 
are that the cold population forms in situ, that a fast 
damping (0.4 Myr as opposed to 3 Myr) or precession 
rate is required, and that resonances, rather than being 
the mechanism for creating the cold population by over- 
lapping to create a chaotic sea, cannot overlie the cold 
classical region while Neptune's eccentricity is high. 

Another key finding is that the "chaotic sea" that may 
have existed in the classica l region during Neptune's wild 
days (jLevison et al.1 120081 ) would not have been caused 
solely by Neptune's high eccentricity but by oscillations 
in Neptune's semimajor axis due to its near-resonant in- 
teractions with Uranus. Thus the exact dynamical con- 
figuration of Uranus and Neptune controls the extent and 
existence of a chaotic zone. Configurations of these two 
planets in which their interactions are especially strong 
might rule out the region of parameter space of 28 AU 
< on < 29 AU during Neptune's high eccentricity period, 
which corresponds to the classical objects being sand- 
wiched between the 9:5 and the 2:1. A detailed investiga- 
tion of the effects of the interactions between Uranus and 
Neptune will likely provide additional constraints on the 
dynamical history of the solar system. The conditions 
for the chaotic sea not to disrupt the cold classicals may 
rule out additional parameter space but will not open up 
more. 

There is a major outstanding problem with forming 
the cold classicals in in situ: unsettlingly, the current 
surface density of cold classicals is thought to be too low 
for the in s itu formation of the 100 — 1000 km object s 
we observe (jStern fc Colwelllll997Hkenvon fc Luulll998D . 
One potential explanation is that the population has lost 
substantial amount of mass to collisions and subsequent 
removal by radiation forces. Another potential resolution 
is that, given that the physics of planetesimal formation 
is currently poorly understood, it may be possible to form 
such large objects at such low surface densities (if some 
major physical process is missing from our understanding 
of planetesimal formation) . Finally, as discussed in Sec- 
tion [531 the cold classical region may have been depleted 
by scattering before the period of hot classical delivery. 
Though none of these potential solutions have been val- 
idated, in situ formation of the cold classicals remains 
viable due to their distinct physical properties. 

Our results are intended to provide constraints for ex- 
tensive numerical integrations that include all the giant 
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planets, have tens of thousands of test particles, and last 
for the age of the solar system. We established conser- 
vative criteria in order to confidently rule out regions of 
parameter space; the remaining regions are potentially 
viable but may be ruled out by additional constraints, 
including those that depend on the details of the config- 
uration of giant planets. We have focused on the classical 
KBOs in this paper and have not tried explicitly to match 
the distribution of resonant KBOs. We expect resonant 
objects to be produced, within the parameter space we 
const rained, by a com bination of migratio n (jMalhotral 
1995), chaotic capture (jLevison et al.ll2008l ). and, a new 
mechanism identified here as being important, fast secu- 
lar evolution to low eccentricities of hot objects delivered 
near resonance followed by capture. Detailed matching 
of the resonant population is beyond the scope of this 
paper but will likely tighten our constraints. 

The scattered disk population may also provide ad- 
ditional constraints on Neptune's dynamical history. 
Within the generalized model we have explored, these 
objects originate in the inner disk and are scattered out 
to beyond 48 AU, or within 48 AU at higher eccen- 
tricities and inclinations than the classical objects. If 
the scattered and hot objects have the same origin, any 
model must correctly produce their relative number. The 
scattered disk also contains a number of objects beyond 
50 AU occupying high-o rder mean motion resonances. 
ILvkawka &: Mukal (|2007D found that in order to produce 
these resonant objects, Neptune must undergo migration 
after the Kuiper Belt has been pre-excited out to 50 AU. 
Our constraints may be consistent with this requirement, 
since we find that Neptune's eccentricity should damp be- 
fore it migrates to its current location and that, in one al- 
lowed region, unexcited objects are sandwiched between 
objects excited by fast secular evolution, extending out 



to the edge of the 2:1 resonance. We leave detailed ex- 
plorations of constraints from the population of resonant 
objects in the scattered disk for future work. 

The fast damping of Neptune's eccentricity would im- 
ply frequent planetesimal scatterings and thus a high sur- 
face density of planetesimals in the vicinity of its orbit. 
Fast precession would imply strong interactions with the 
other giant planets (or a high surface density disk). Ei- 
ther way, the dynamical histories of Neptune that pro- 
duce the hot and cold KBOs are very different from the 
peaceful disk formation that was the paradigm until 15 
years ago. The most viable regions of parameter space 
we identified imply that Neptune underwent a period of 
high eccentricity but, mercifully for the spared planetes- 
imals that are today's cold classicals, Neptune's "wild 
days" were over soon. 
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APPENDIX 

STATISTICAL SIGNIFICANCE OF THE HOT AND COLD CLASSICAL ECCENTRICITY DISTRIBUTIONS 

In addition to the qualitative assessment performed in Section [21 we also conduct statistical tests of the significance 
of the confinement of the cold population. We perform the one-dimensional Kolmogorov-Smirnov (KS) test with the 
null hypothesis that the observed distribution is consistent with being drawn uniformly in a and e and then filtered 
by the survival map (jLykawka & Mukai 2005). We perform the test separately in two regions, 42 < a < 44 AU and 
44 < a < 45 AU, because the survival map is different in these two regions. Note that we use 42 AU instead of 
42.5 AU in order to increase the sample size. We created the survival map distribution from an initial distribution 
uniformly spaced in e < 0.3 with the "stability map filter" method described in Section [23] for three assumed survival 
rates: bottom of the range, middle of the range, and top of the range. We increased the number of objects in the 
survival map distribution until the results converged, which meant we had to use 10,000 objects post stability map 
filtering. First, for the observed distribution, we compared to the observed "likely cold" objects with i < 2°. The 
resulting probabilities that the observed distribution is consistent with a population shaped only by long-term stability 
are summarized in Table [TJ 

TABLE 1 

Probability from KS t est comparing observed cold objects with i < 2° from Minor Planet Center (IGladman et al.] 2008; 
IVolk & Malhotra 20TJJ) to "surv ival map" distribution. " Low," "mid," and "high" refer to the bottom, middle, and top of 

THE 10% SURVIVAL RANGE USED BY LYKAWKA <fc MlJKAI (2005)). FOR EXAMPLE, FOR THE SURVIVAL RANGE OF 50 — 60%, LOW, MIDDLE. 
AND HIGH WOULD INDICATE THAT 50%, 55%, AND 60% OF PARTICLES IN THAT (a, e) CELL SURVIVE. SEE SECTION 12.31 AND APPENDIX lAl FOR 

DETAILS. 

sample size low mid high 
42 AU < a < 44 AU 38 10~ 7 10~ 8 10~ 8 

44 AU < a < 45 AU 13 0.004 0.002 0.001 

Then we created five alternative samples of cold objects. Instead of choosing objects with i < 2° for inclusion in the 
cold sample, we selected a uniform random number between and 1 for e ach observed object . If the number was less 
than the probability that the object is cold (based on the distribution of iGulbis et al.l ((2010)), we included it in the 
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sample size 


low 


mid 


high 


60 


lO" 8 


io- 10 


io- 11 


64 


io- 10 


io- 11 


io- 13 


69 


10 -io 


io- 12 


io- 13 


67 


10~ 9 


io- 11 


icr 12 


69 


HT 11 


io- 12 


1Q -14 


31 


IO" 5 


10~ 5 


IO" 6 


33 


HT 5 


10~ 5 


IO" 6 


30 


KT 4 


10~ 5 


10~ 6 


32 


IO" 4 


IO" 5 


IO" 5 


29 


IO" 5 


IO" 6 


IO" 6 



sample. Using the iBrownl (|2001l) distribution of inclinations instead of the iGulbis et al.l ()201ClO did not significantly 
affect our results. The resulting probabilities are summarized in Table [2j 

TABLE 2 

Probability from KS test comparing probabilitistically-selected observed cold objects to "survival map" distribution. 
42 AU < a < 44 AU 



44 AU < a < 45 AU 



Based on these results, the orbital distribution of observed objects is not consistent with our null hypothesis. The 
confinement of cold objects to low eccentricities is formally statistically significant. 

Then we repeated the tests using only the CFEPS objects. Instead of using the survival map distribution, we took 
the survival map distribution and applied the CFEPS Survey Simulator. We assumed either H magnitudes uniformly 
distributed between 6 and 8 or randomly selected from the observed classical objects, and the results were insensitive 
to this choice. First we used the "likely cold" objects with i < 2°. The resulting probabilities are summarized in Table 

13 

TABLE 3 

Probability from KS test comparing observed CFEPS cold objects to survey-simulated "survival map" distribution. 

sample size low mid high 
42 AU < a < 44 AU 12 0.13 0.06 0.04 

44 AU < a < 45 AU 11 0.09 0.06 0.03 

Then we created five alternative samples of cold objects, as described above for Table 2. The resulting probabilities 
are summarized in Table 3J 

TABLE 4 

Probability from KS test comparing observed CFEPS cold objects to survey-simulated "survival map" distribution. 
42 AU < a < 44 AU 



44 AU < a < 45 AU 



sample size 


low 


mid 


high 


20 


0.01 


0.003 


0.002 


25 


0.001 


0.0004 


0.0002 


20 


0.01 


0.003 


0.002 


20 


0.01 


0.003 


0.002 


22 


0.001 


0.0003 


0.0001 


16 


0.03 


0.01 


0.005 


20 


0.009 


0.004 


0.001 


20 


0.009 


0.004 


0.001 


19 


0.01 


0.005 


0.002 


18 


0.02 


0.007 


0.003 



These probabilities using the small sample of CFEPS objects are low, supporting our conclusion that the cold objects 
are confined to low eccentricities, but this result is statistically marginal. However, we would not expect the observed 
objects in the full MPC sample to preferentially have lower eccentricities (indeed, in Figure |H the Survey-Simulated 
survival map distribution follows the survival map closely), so, given our results for the full MPC sample, we expect 
the significance to increase as the CFEPS sample size becomes larger. 

An alternative statistical test is the Anderson-Darling test, which is more sensitive to the tail of the distribution. 
However, even though we took measures to avoid contamination, there are likely to be "contaminating" objects in the 
observed cold objects that are actually hot. Therefore we do not necessarily want to give the outliers higher weight, 
so we judge that the KS test is more robust for this purpose. Using the Anderson-Darling test, we obtained similar 
results in the cases with large sample sizes and somewhat higher probabilities in the low sample size cases. 

Proper Elements 

We now c onsider the free, or prop er , elements of t h e obs e rved KBOs, which have b een computed for a subset 
of KBOs by IKnezevic fc Milanil (|2000f ): IKnezevic etall (|2002fl ; IKnezevic & Milanil (|2003t ). The free elements precess 
about the forced values, which are set by the current configuration of the giant planets, and thus provide a better 
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window to the history of the solar system than the instantaneous orbital elements. In Figure 1251 we plot the prope r 
element eccentricities and inclinations of observed KBOs on top of the survival maps of Lykawka & Mukai (2005), 
which are formulated in terms of instantaneous eccentricity and inclination. Qualitatively, we see the same features 
as in Figure |3]and Figure |U the cold classicals (red squares) are confined below e < 0.1 in the region from 42.5 to 45 
AU. Throughout the region, the hot objects (blue triangles) occupy the upper portion of the stability region. 




42 44 46 48 42 44 46 48 

a (AU) a (AU) 

10 20 30 40 50 60 70 80 90 100 
Survival Rate (Lykawka and Mukai 2005) 

Fig. 25. — Plotted over the survival maps of Lvkawka & Mukai (2005) are the proper eccentricity (left) and proper inclination (right) 
distributions of the observed classical KBOs. Note that the survival map uses instantaneous orbital elements. The red squares are objects 
with i < 2° and are thus very likely cold classicals. The blue triangles have i > 6° and are thus very likely hot classicals. The membership 
of any given purple circle (2° < i < 6°) is ambiguous. These inclinations are now defined using the proper elements. The light red squares, 
light blue triangles, and light purple circles, respectively, are objects for which proper elements have not been comp uted, and thus we plot 
their instan taneous elements . Classical objects are taken from the Minor Planet Cente r Data b ase and classifie d by Gladrnan ct al. (2008) 
and Volk & Malhotra (2011). Proper elements were computed by Knczcvic & Milani (2000); Knczcvic ct al. (2002); Knczcvic & Milani 
(2003). The yellow line indicates our conservative criterion for preserving the cold classicals. 



We repeat the same statistical tests (Table [5] and [6]) as above for the subset of observed KBOs that have computed 
proper elements. We use the proper inclinations to classify the objects with i < 2 a s cold. For the second t est, for 
which the sample is probabilistically-selected, we use the inclination distribution from iVolk fc Malhotral (pOll . which 
considers the inclinations with respect to the invariable plane. The results are consistent with the those from the 
instantaneous orbital elements above. We caution that we are comparing the proper eccentricities of observed objects 
to a stability map of instantaneous eccentricities. This comparison should be repeated when a stability map formulated 
in terms of proper elements becomes available. 

TABLE 5 

Probability from KS test comparing observed proper eccentrkitie s cold objects with i < 2° from Minor Planet Center 
(IGladman et AL.Il2008t IVolk Malhot ra 2011) to "survival map" distribution. 

sample size low mid high 
42 < a < 44 AU 35 10~ 8 10" 9 10" 9 

44 < a < 45 AU 25 0.002 0.0006 0.0003 



DERIVATION OF SECULAR THEORY 

In Section IB. 11 we derive additional factors that we include in the secular theory (Section 14. 2p . In Section IB. 21 we 
consider the secular forcing due to other planets besides Neptune and demonstrate that, for KBOs, the secular forcing 
due to all four planets reduces to the forcing by a precessing Neptune. 
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TABLE 6 

Probability from KS test comparing probabilistically-selected observed cold objects (using proper elements) to "survival 

map" distribution. 



sample size 
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mid 


high 
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io- 11 


IO" 12 


55 


io- 11 


IO" 13 


IO" 13 


51 


io- 12 


IO" 13 


1Q -14 


19 


io- 10 


io- 11 


io- 12 


40 


IO -5 


10" 6 


IO' 6 


37 


10" 4 


10" 4 


IO -5 


39 


io- 4 


10~ 5 


IO" 5 


37 


io~ 4 


10~ 4 


IO' 5 


36 


io~ 4 


io -5 


10~ 6 



Derivation of additional terms 

In Section 14.2.21 we r elegated to the appendix the derivation of several additional terms in the modified secular 
theory. In Section [B. 1.11 we derive the factors proportio nal to e^j that ap pear in the extra factors used in Equation ||5j 
and ([5]) as coefficients to e N . In Section IB. 1.21 we follow Malhotra et al.] (|1989t ) to derive resonance correction terms. 

High order eccentricity terms 

The basic secular theory includes only the lowest order eccentricity terms. However, when Neptune's eccentricity 
is high, terms containing e^ are no longer negligible. Therefore <?kbo and e~f orce d must be modified. Here we define 
the extra ter ms and factors used in E quation ([5]), which come from additional terms in the disturbing function (see 
Chapter 7 of iMurrav fc Dermottll2000l for a standard derivation). The disturbing function has the form, up to second 
order, of: 

i? = nV^[ e 2 (/ 2 + / 5 4 + / 6e 2 ) 

TO 

O 

+ cos(tn - n7 N )ee N (/i + A^n + / 12e2 ) + cos(2(tu - w N ))e 2 e^f 17 ] 

(Bl) 

Fully incorporating all e and eN to second order would modify the functional form of the secular theory. However, if 
we treat eN as a constant, we can modify the f% term in the secular forcing frequency ^kbo to /a + /se 2 ^ (Equation 
[5]) and the /io term in the forced eccentricity to f w + /ne N (Equation [8]). Because the form of the secular evolution 
Equation ([2]) is derived by differentiating R with respect to the particle's h and k, treating eN as a constant does not 
modify the form o f the secular evolution equ ations but simply adds extra correction factors. The / factors are defined 
in Appendix B of IMurrav fc Dermottl ((2000). The success and necessity of these extra terms is illustrated in Figure 

m 
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Fig. 26. — Evolution of several example particles illustrating the necessity of the additional e 2 ^ terms. The colored lines, each representing 
one of three particles, are the output of a numerical integration in which Neptune has ejj = 0.3. Each particle has a different linestyle. 
The black curve is the analytical model with higher order eN terms incorporated (Equation [4} while the gray curve is the analytical model 
that neglects these higher order terms (Equation |5J. For each particle, the black curve matches better than the gray curve. When the 
particle itself has a high eccentricity, depending on its initial phase, the frequency of the analytical model matches well (purple) or is off 
by a factor of up to two (cyan worst case scenario). The discrepancy between the analytical model and the numerical integration output 
at some phases for high eccentricity particles is the result of the approximation which is only to lowest order in the particle's eccentricity. 
Fortunately, the most favorable case for delivering the hot classicals — an initial periapse that results in the particle's eccentricity sharply 
decreasing — is that for which our analytical model performs best. 
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Modification near resonances 

Proximity to resonance modifies the secular fre quencies. Following IMalhotra et al.l (|1989l ). we define the resonance- 
correction terms C and e in Equations (5) and (6). IMalhotra et all (|1989f) developed equations for a pair of moons near 
first order resonance. Instead of a pair of moons, we treat a massless KBO and Neptune, and we include resonances 
above the first order. 

Resonances add additional important terms to the disturbing function, R. The factor C in Equation (6) is propor- 
tional to the coefficient of the direct part of the disturbing function, with argument (j + x)X — j'An — xm^, where 
x is the order of the resonance. Since the coefficients of this part of the disturbing function R are proportional to 
efj, and since h oc ^ and k oc a factor of xe^f 1 comes in (Equation 6), which was not explicitly included in 
IMalhotra et al.l (|1989f ) because they treated only the x = 1 case. The factor C for each order x i s tabulated in Table 
[71 Th ese coefficients were taken from the expansion of the disturbing function in Appendix B of iMurrav &: Dermotti 



TABLE 7 

Coefficients for Equation (6). 

x C 

1 !(-2Ci+*)-a£)6g+«>(«) 

2 |((-5(j + x) + 4(j + xf) + (-2 + 4(j + x))a£ 

3 ^((-26(j + x) + 30(j + xf - 8(j + x) 3 ) + (-9 + 27(j +x) - 12(j + xf)a£ 

4 ^((-2060' + x) + 283(i + xf - 120(j + xf + 16(j + xf) + (-64 + 236(j + x) 
-I68(j + xf + 32(j + xf)a£ 

+(48 - 78(j + x) + 24(j + xfWj^ + (-12 + 8(j + x))a 3 

+« 4 ^)^ / r ) («) 



The factor e in Equation ([3]) depends on the proximity to resonance. Extending IMalhotra et all (|1989f) to resonances 
of arbitrary order, we obtain: 

w = jn N - (j + x)n(l + -^-(1 + a ^-)bf) ) 

m da 1 

O 

t=l — U + *) 2 2 7— 7^2 

2 m (u> n) 

O 

(B2) 

where njy is the mean motion of Neptune and n is the mean motion of the particle. 

The n : 1 resonances have indirect terms not explicitly included in Malhotra et al.l ((1989) (R. Malhotra, private 
communication). However, the only relevant n : 1 resonance in the region of the Kuiper Belt we are studying is the 
2 : 1 resonance and its indirect terms result in expressions that, when incorporated above, are directly proportional to 
the mass of the KBO and thus assumed to be negligible. 

Effects of other planets 

In the case of multipl e planets, the forced eccentricity of a small body on an external orbit is given by 
(M urrav fc Dermottl[200a i: 

N 

ho = ~ -j— ! — sm(git + Pi) 

N 

k = - V — - — cos(g;i + Pi) 
i=i A ~ 9t 



(B3) 



where 



3=1 
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1 TOo (2) , 

o 

N 1 

i=i c 

(B4) 

where the particle's forced eccentricity ef orcc d = V^o + ^o an d gi and ey, are the eigenfrequencies and eigenvector 
components of the planetary system. We compare Equation (|B3|) with Equation (7). For KBOs in the classical region, 
Neptune's orbits dominates A, the precession rate of the particle's free eccentricity. The precession rate of a particle's 
free eccentricity due to Neptune alone agrees with the four-planet case to within 30%. Thus A ps <?kbo Neptune 
dominates the Vi term, and thus the four-planet secular theory reduces to the single-planet secular theory with an 
extra <ft = t^n term for Neptune's precession. Therefore the four planet Equation (|B3|) reduces to Equation (7). 
The quantity z/j corresponds to 3KBO e N(^)> 9i corresponds to zujsi, and A corresponds to <?kbo- The extra sin term in 
Equation (7) is an empirical factor to acc ount for eccentri c ity da mping and is discussed in detail in the main text. This 
conclusion is consistent with the result of iChiang fc Choi ([2008) that the current forced eccentricities of the KBOs are 
largely determined by Neptune's orbit. 
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